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ABSTRACT 


Refraction  data  along  the  structural  axis  of  Vancouver 
Island  were  recorded  as  part  of  the  1980  Vancouver  Island 
Seismic  Project.  The  refraction  line  extends  about  330  km 
from  the  southern  to  the  northern  end  of  the  island.  These 
partially  reversed  data  were  interpreted  by  the  analysis  of 
first  arrivals  in  terms  of  a  horizontal,  plane-layered  model 
and  also  using  a  WHB  inversion  of  travel  times.  Synthetic 
seismograms  were  also  computed  in  an  attempt  to  match 
possible  Moho  reflections. 

Interpretation  of  these  data  suggests  an  increasing 
velocity  from  5.30  km/s  at  the  surface  to  6.5  km/s  at  5  km 
depth,  and  a  more  gradual  increase  to  7.0  km/s  at  17  km 
depth.  Discontinuity  of  the  travel-time  curve  on  the  Fuca 
profiles  has  been  interpreted  as  evidence  of  a  low  velocity 
zone  in  the  lower  crust.  Although  probable  Pn  arrivals  from 
the  Moho  were  observed,  its  depth  and  shape  remain 
uncertain.  However,  an  interpretation  that  fits  the  observed 
refraction  data  from  the  crust  and  reflection  data  from  the 
Moho  by  synthetic  modelling,  has  a  Moho  depth  of  44  km  and  a 
mantle  velocity  of  8.0  km/s.  Lateral  variation  of 
thicknesses  is  significant  beneath  the  central  part  of 
Vancouver  Island.  The  analysis  of  first  arrivals  based  on 
the  limited  refraction  data  makes  the  interpretation 
somewhat  arbitrary,  thus  the  proposed  model  is  only 
speculative  at  the  present  time. 
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1.  INTRODUCTION 


The  study  of  transmission  of  elastic  waves  through  the 
earth  has  been  one  of  intensive  research  in  geophysics  over 
the  past  few  decades.  Seismology  provides  the  tools  for 
studying  the  structure  and  composition  of  the  earth’s 
interior.  Seismic  evidence  has  shown  that  the  crust  of  the 
earth  is  about  50  km  thick  under  continents  and  5  km  under 
the  oceans;  a  distinct  worldwide  discontinuity  called  the 
'Moho'  (the  Mohorovicic  discontinuity),  which  has  the 
characteristics  of  strong  seismic  velocity  variations, 
divides  the  crust  from  the  upper  mantle. 

In  Canada,  deep  crustal  studies  started  in  1947 
(Hodgson,  1947).  Recently  ten  groups  from  Universities  and 
Government  have  carried  on  a  long-term  seismic  program  to 
study  the  crust  and  upper  mantle  over  different  regions. 
Some  examples  of  the  studies  that  were  done  included: 
Cumming  and  Kanasewich  (1966);  Rankin,  Ravindra  and 
Zwicker ( 1 969  )  ;  Mereu  and  Hunter (  1 969 ) ;  Overton  (  1970); 
Barr(1971);  Hales  and  Nation  (1973);  Berry  and 
Forsyth ( 1 975 ) ;  Cumming,  Clowes,  and  Ellis( 1979) . 

A  review  of  the  published  interpretations  leads  to  a 
general  crustal  structure  in  Canada  ( Berry , 1 973 ) : 

1.  Pn  velocities  ( Compress i onal  wave  velocities  under  the 
Moho  discontinuity)  vary  from  a  value  possibly  as  low  as 
7.7  km/sec  under  Vancouver  Island  to  8.6  km/sec  and 
higher  in  the  extreme  eastern  part  of  the  shield  and 
some  parts  of  Atlantic  coast. 
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2.  Large  areas  of  Canada  have  a  crustal  thickness  of  30-40 
km, with  Vancouver  Island,  the  southwestern  Prairies,  the 
Lake  Superior  basin  and  parts  of  the  eastern  shield  of 
Quebec  being  thicker. 

3.  The  Riel  discontinuity,  a  deep  intra-crustal  reflector 
and  sometimes  a  refractor,  is  widely  reported  in  the 
Prairies,  and  Arctic  polar  continental  shelf.  It  is  not 
seen  in  Hudson  Bay,  nor  under  Lake  Superior,  nor  in 
Maritime  regions. 

This  investigation  is  part  of  large  scale  seismic 
program  (VI SP  80)  which  consists  of  both  seismic  refraction 
and  reflection  experiments  in  the  region  of  Vancouver 
Island.  The  overall  objective  of  this  seismic  program  is  to 
explore  the  nature  of  the  crust  and  upper  mantle  in  the 
transition  zone  from  the  Pacific  Ocean  through  Vancouver 
Island  to  the  mainland  of  British  Columbia. 

The  present  study  has  the  following  purposes: 

1.  Develop  crustal  velocity-depth  functions  for  refraction 
data  collected  from  Vancouver  Island. 

2.  Utilize  these  velocity-depth  functions  for 

interpretation  of  Vancouver  Island  crustal  structure, 
and  provide  a  control  point  on  the  cross  refraction 
profiles  which  will  be  interpreted  by  others. 


2.  BACKGROUND  GEOLOGY  AND  GEOPHYSICS 


2.1  Generalized  Geology  of  Vancouver  Island 

Vancouver  Island  is  located  in  the  Insular  Belt  of 
western  British  Columbia.  The  Insular  Belt  is  a  series  of 
complexes  of  largely  igneous  (volcanic)  origin,  believed  to 
have  formed  as  part  of  an  island  arc  system  in  early 
Mesozoic  time.  They  are  probably  allochthonous  (Monger  and 
Price,  1979).  The  events  of  mountain-building  in  Mesozoic 
time  had  a  significant  impact  on  the  formation  of  Vancouver 
Island.  These  events  involved  widespread  folding,  faulting, 
and  the  intrusion  of  batholiths  in  the  Insular  Belt  itself, 
the  Coast  Mountains,  and  the  Intermontane  Belt. 

Vulcanism  on  Vancouver  Island,  in  the  form  of  basaltic 
lava-flows  was  extensive  during  the  Mesozoic  age.  As  a 
result,  thick  sequences  of  volcanic  rocks  of  oceanic  crustal 
and  island  arc  type  occur  throughout  the  island.  Erosion  and 
differential  uplift  continued  to  the  end  of  the  Tertiary 
period.  The  complex  history  of  Vancouver  Island  has  resulted 
in  a  great  variety  of  volcanic,  plutonic,  sedimentary,  and 
metamorphic  rocks  which  result  in  a  high  variability  of  near 
surface  lithologies  (Muller,  1974,  1977).  However,  these 
rocks  are  in  general  of  comparatively  high  density  and 
highly  indurated  and  as  a  consequence,  near  surface  seismic 
velocities  are  relatively  high,  compared  for  example  to  the 
plains  of  Western  Canada. 
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2.2  Tectonic  Setting  of  Vancouver  Island 

The  essential  features  of  plate  tectonic  environment  in 
coastal  British  Columbia  are  dominated  by  the  relative 
motion  of  three  main  lithospheric  plates:  the  large  Pacific 
and  America  plates,  and  the  Juan  de  Fuca  plate.  A  small 
northern  part  of  the  Fuca  plate,  which  has  been  named  the 
Explorer  plate,  has  been  shown  to  be  moving  independently. 
Fig. 2.1  shows  the  tectonic  map  in  this  area.  Two  distinct 
features  are  present  along  the  continental  margin:  (Keen  and 
Hyndman , 1979) 

1.  From  the  southern  end  of  the  Vancouver  Island  to  just 
south  of  Queen  Charlottle  Islands  there  is  convergence 
with  subduction  of  ocean  crust  beneath  the  continent. 

2.  From  the  Queen  Charlottle  Islands  to  southern  Alaska 
there  is  a  transform  strike-slip  fault  separating  the 
Pacific  and  America  plates. 

The  motion  between  the  Juan  de  Fuca  and  Explorer  plates 
occurs  across  the  Nootka  transform  fault  zone,  which  extends 
northeasterly  from  the  north  end  of  the  Fuca  ridge  to  the 
continental  shelf  off  north-central  Vancouver  Island.  The 
active  portion  of  the  fault  zone,  about  20  km  wide,  has 
produced  extensive  disturbance  in  the  0.5  to  1  km  of 
overlying  sediments.  Hyndman,  Riddihough,  and  Herzer  (1979) 
showed  that  the  magnetic  anomaly  pattern  required  the  fault 
to  have  a  strike-slip  motion  of  about  3  cm/yr,  with 
subduction  at  the  margin  being  more  rapid  to  the  south  than 
to  the  north  of  the  fault  zone. 
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Figure  2.1  Tectonic  map  of  western  Canada  showing  the  main 
lithosphere  plate  boundaries  and  relative  plate  motions.  The 
dotted  line  is  the  estimated  northern  edge  of  the  Juan  de 
Fuca  plate.  PRfz=Paul  Revere  fracture  zone;  Sfz=Sovanco 
fracture  zone;  PP=Pacific  plate;  EP=Explorer  plate;  JP=Juan 
de  Fuca  plate;  AP=America  plate.  (modified  from  Keen  and 
Hyndman , 1979) 


KM 

Figure  2.2  Structural  section  along  line  CD  in  Fig.  2.1. 
Gravity  profiles  are  free  air  over  sea  and  Bouguer  over 
land;  sol id=observed ,  da shed=computed .  Density  in  g.cm' 3 . 
(after  Riddihough,  1979) 
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Most  subduction  zones  in  other  regions  are 
characterized  by  a  deep  trench,  but  there  is  no  such  deep 
trench  in  this  region.  With  the  evidence  of  a  lack  of  a  deep 
margin  trench  and  deep  earthquakes  in  this  zone,  Riddihough 
(1978)  suggested  that  the  subduction  might  be  a  very  recent 
geological  feature,  and  the  absence  of  a  deep  trench  was 
probably  due  to  high  sedimentation  rates  and  the  damming 
effect  of  the  Juan  de  Fuca  Ridge  lying  only  a  few  100  km 
off  shore . 

Compressional  folding  and  uplift  of  Tertiary  sediments 
are  also  consistently  observed  along  the  coast  regions,  and 
in  detail  a  model  of  imbricate  thrusting  along  the 
continental  margin  fits  the  pattern  of  deformation  closely. 
Detailed  experiments  along  the  continental  margin  attempting 
to  locate  microseismicity,  have  proved  negative.  However, 
there  is  a  strong  concentration  of  seismicity  of  magnitudes 
up  to  6  on  the  Richter  scale  in  the  southern  Georgia 
Strait-Puget  Sound  area.  The  larger  events  generally  are 
deeper  than  normal,  most  at  about  50  km.  Crosson  (1972), 
Roger(1979),  Keen  and  Hyndman ( 1 979 )  concluded,  from  the 
earthquake  study,  that  there  were  two  distinct  groups  of 
earthquakes:  a  shallow  group  and  a  second  deeper  group.  An 
inclined  zone,  dipping  at  a  shallow  angle  and  along  which 
oceanic  lithosphere  is  supposedly  subducting,  was  postulated 
to  divide  the  two  groups  of  earthquakes.  Thus  the 
concentration  of  earthquakes  in  the  Georgia  Strait  region 
indicates  that  oceanic  lithosphere  moves  under  the 
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continental  margin  at  a  shallow  angle  until  the  region  of 
Georgia  Strait  and  Puget  Sound,  then  dips  steeply  about  10 
to  15  degrees  down  to  the  mantle. 

Heat  flow  and  gravity  patterns  provide  additional 
evidence  of  subduction  in  this  region.  The  heat-flow  has  the 
characteristic  pattern  of  a  subduction  zone:  a  band  of  low 
heat  flow  extends  from  the  trench  to  the  volcanic  arc,  and  a 
much  higher  than  normal  heat  flow  extends  for  a  considerable 
distance  inland  (McKenzie  and  Sclater  1968;  Hyndman , 1 976 ) . 
The  low  zone  is  produced  by  the  cold  oceanic  lithosphere 
that  sinks  beneath  the  continent;  the  high  heat  flow  farther 
inland  probably  results  from  the  sinking  slab.  The  gravity 
anomaly  exhibits  a  similar  characteristic  pattern  of  two 
pairs  of  low  and  high  bands,  as  shown  in  Fig.  2.2.  The  main 
inland  band  reflects  the  established  axis  of  subduction  of 
the  Juan  de  Fuca  plate;  the  weaker  band  at  the  continental 
margin  appears  in  part  to  be  produced  by  the  effect  of  the 
boundary  between  thin  oceanic  and  thick  continental  crust. 

Riddihough  (1979)  developed  a  detailed  structural  model 
from  gravity  data.  Fig. 2. 2  shows  his  model  across  northern 
Vancouver  Island.  The  model  suggests  that  oceanic 
lithosphere  is  descending,  bending  and  disappearing  near  the 
base  of  the  continental  slope,  and  that  the  crust  under 
Vancouver  Island  is  between  20  to  30  km  depth  and  is 
underlain  by  anomalous  high  density  material  to  produce  the 
observed  gravity  field.  On  the  whole,  all  geophysical 
evidences  strongly  support  the  existence  of  the  subduction 
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zone  in  southwest  British  Columbia. 


2.3  Previous  Seismic  experiments 

Vancouver  Island  has  been  the  subject  of  considerable 
geophysical  investigations.  The  early  structural 
interpretation  was  based  on  the  refraction  work  of  White  and 
Savage  (1965),  and  time-term  analysis  of  the  same  data  by 
Tseng  (1968).  White  and  Savage  suggested  an  unusually  thick 
crustal  section  under  Vancouver  Island  with  the  Moho  lying 
at  a  depth  of  50  km  and  with  a  Pn  velocity  of  7.7  km/s. 
Tseng  obtained  a  similar  crustal  thickness,  but  he 
introduced  an  intermediate  layer  at  a  depth  of  25  km  with  a 
P-wave  (compressional  wave  )  velocity  of  7.1  km/s. 

Berry  and  Forsyth  (1975)  attempted  to  reconcile  their 
data  with  earlier  work.  Their  data  indicated  a  major 
discontinuity  in  the  travel-time  curve  between  the  mainland 
coast  and  Vancouver  Island.  This  major  discontinuity  is 
likely  to  represent  a  zone  of  scattering  for  seismic  energy 
in  that  region.  They  also  found  a  thick  crust  of  about  50  km 
and  an  intermediate  layer  at  about  30  km  under  Vancouver 
Island.  Furthermore,  analysis  of  surface  waves  by 
Wicken(1977)  concluded  that  two  velocity  discontinuities 
existed  at  about  these  depths.  He  indicated  that 
inconsistencies  between  the  model  values  for  Love  and 
Rayleigh  waves  were  suggestive  of  upper  mantle  anisotropy. 
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However,  from  the  analyses  of  P-S  conversions  of 
long-period  teleseismic  P-waves  recorded  at  the  southern  end 
of  Vancouver  Island,  Langston ( 1 98 1 )  suggested  that  there 
were  two  Moho  discontinuities  in  this  region.  The  first  Moho 
was  at  a  depth  of  20  km,  in  which  the  material  at  this  depth 
represented  a  remnant  of  high  density  upper  mantle,  and  the 
second  was  at  a  depth  of  45  km  in  which  the  Moho  was 
interpreted  as  the  subducting  oceanic  crust  of  the  Juan  de 
Fuca  plate.  A  low  velocity  zone  between  20  km  and  45  km  was 
also  postulated  to  divide  the  two  Moho  discontinuities.  If 
the  low  velocity  zone  exists  in  this  region,  it  would  no 
doubt  create  difficulties  in  interpretation  of  refraction 
work  . 

Riddihough( 1 979 )  provided  a  summary  of  all  seismic  and 
gravity  interpretations  beneath  Vancouver  Island.  A  feature 
of  these  interpretations  is  the  presence  of  two 
discontinuities,  one  near  30  km  depth  and  one  at  between  40 
to  50  km  depth.  The  crustal  thickness  determined  by  seismic 
methods  is  significantly  different  from  that  deduced  by 
gravity  data.  To  reconcile  the  conflicts  between  two 
interpretations,  Riddihough  suggested  that  the  subducting 
oceanic  crust  of  the  Juan  de  Fuca  plate  beneath  Vancouver 
Island  was  between  30  to  50  km.  The  material  of  this 
subduction  zone  has  a  special  feature  of  having  high  density 
and  low  P-wave  velocity.  Fig. 2. 3  shows  the  P-wave 
velocity-depth  interpretations  proposed  in  a  number  of  these 
studies.  The  wide  range  of  the  interpretations  in  Fig.  2.3 
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Figure  2.3  P-wave  velocity  profiles.  These  are  the  results 
of  some  previous  studies  of  structure  beneath  Vancouver 
Island  and  adjacent  regions  in  the  United  States:  solid 
line=  Langston,  dashed  line=  Berry  &  Forsyth,  stars=  White  & 
Savage,  (after  McMechan  and  Spence,  1982) 


suggests  that  either  the  data  are  ambiguous  or  there  are 
significant  lateral  variations  in  structure  beneath 
Vancouver  Island. 

The  seismic  data  existing  prior  to  this  project  do  not 
define  clearly  the  structural  characteristics  of  the  crust 
and  upper  mantle  in  this  region.  The  present  study  is  part 
of  a  continuing  investigation  aimed  at  clarifying  the 
structure  under  Vancouver  Island. 


3.  DATA  ACQUISITION 


3.1  Location  of  Refraction  Lines 

In  August  1980  the  CO-CRUST  group  conducted  a  field 
program  on  Vancouver  Island,  and  included  workers  from  the 
Universities  of  British  Columbia,  Alberta,  Saskatchewan, 
Manitoba  and  Western  Ontario,  and  the  Department  of  Energy, 
Mines  and  Resources  of  Ottawa.  The  initial  phase  of  the 
experiment  consisted  of  refraction  measurements  undertaken 
in  the  period  August  9  to  21  followed  by  a  reflection 
experiment  from  August  22  to  26. 

The  data  for  the  present  study  consist  of  two 
refraction  lines,  NA  and  AF  in  Fig. 3.1.  The  distance  between 
the  two  outermost  locations,  N  and  F  is  approximately  330 
km.  The  centre  of  the  refraction  line  is  located  at  A. 
Initially  38  seismic  systems  were  distributed  along  the  150 
km  section  N-A.  Shots  A1,  N1,  and  FI  were  then  detonated  at 
A , N , F  respectively.  The  seismometers  were  relocated  along 
180  km  of  the  line  A-F  for  the  observation  of  shots  F2,  A2, 
and  N2  at  F,  A,  N  respectively. 

Each  of  these  six  profiles  is  referred  to  by  a  name 
consisting  of  a  shot  point  designation  (N,A,  or  F)  followed 
by  a  number  ’1'  or  ’2’.  '1T  refers  to  the  northern  profile; 
’2'  to  the  southern  profile.  For  example,  profile  N1  refers 
to  a  shot  at  N  recorded  along  the  northern  line  NA.  The 
charge  size  for  six  shots  varied  from  900  to  1800  kg.  The 
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Figure  3.1  VISP  80  Location  map.  The  data  are  from  three 
shot  points  N,  A,  and  F.  Portable  seismometers  were  placed 
along  the'solid  lines  at  an  average  spacing  of  about  5 
km. (from  McMechan  and  Spence,  1982) 
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distant  shots  were  1800  kg;  the  nearer  shots,  900  kg.  For 
details  of  shot  point  parameters,  the  reader  is  referred  to 
Ellis  and  Clowes  (1981, p. 12). 


3.2  Data  Characteristics  on  Tape 

All  refraction  data  were  recorded  on  magnetic  tape  by 
eight  different  types  of  recording  systems,  five  digital  and 
three  analog.  Each  system  has  its  own  characteristic  gain 
and  frequency  response  to  record  seismic  signals  over  a 
passband  from  1  to  20  Hertz.  Ellis  and  Clowes ( 1 98 1 , p . 20 ) 
discussed  the  details  of  each  system.  Each  seismometer 
recorded  at  least  60  seconds  of  signals  after  the  first 
arrival  of  energy.  WWVB  signals  were  also  recorded  at  both 
the  recording  site  and  shotpoint  to  maintain  time  control. 

All  field  data  were  redigitized  at  a  rate  of  60  Hz  with 
variations  in  analog  tape  speeds  taken  into  consideration. 
Each  record  contains  7200  samples,  and  the  first  three  data 
words  provide  information  concerning  receiver  number,  shot 
number,  and  orientation  of  the  seismometer.  The  data  were 
written  in  format  (5E16.6)  with  a  block  length  of  1200 
samples  (4800  bytes).  Thus  each  event  consists  of  6  blocks 
of  data  on  the  magnetic  tape.  To  convert  all  the  data  to  a 
common  ground  velocity,  each  recording  has  been  modified  by 
a  gain  factor  which  reflects  the  various  sensitivies  of  the 
seismometers  and  the  gain  setting  of  the  recording 
ampl i f ier s . 
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3.3  Instrumentation  of  the  University  of  Alberta 

In  1980  the  seismic  array  recording  system  of 
University  of  Alberta  had  twelve  channels  to  record  seismic 
data  and  two  additional  channels  to  record  shot  instant  from 
a  radio  receiver  and  absolute  time  from  WWVB .  All  data  were 
recorded  on  a  Peripheral  Equipment  Corporation  model  7830-9 
digital  tape  recorded  at  a  tape  velocity  of  6.25  inches  per 
second,  yielding  a  transfer  rate  of  5000  bytes  per  second. 
(Alsopp,  Burke,  and  Cumming  1972) 

The  impulse  response  of  the  recording  instrument  and 
seismometer  can  be  derived  from  an  S-plane  representation 
(Dickins , 1 973 ) .  The  constants  in  the  S-plane  representation 
of  both  recording  system  and  seismometer  may  be  found  by 
examining  the  specification  of  their  circuits  and  their 
components.  Then,  these  nominal  values  may  be  adjusted 
slightly  so  that  the  amplitude  and  phase  spectra  calculated 
from  the  specifications  match  the  amplitude  and  phase 
spectra  determined  experimentally.  The  constants  in  the 
S-plane  representation  were  derived  by  G.L.  Cumming 
according  to  this  procedure  (in  Appendix).  Fig.  3.2  and  3.3 
show  the  amplitude  and  phase  spectra  of  the  amplifier  system 
and  seismometer  which  were  derived  from  the  S-plane 
representation . 

The  amplifier  has  a  flat  response  from  0.1  Hz  to  50  Hz, 
and  the  response  of  the  seismometer  increases  as  the 
frequency  increases.  Alsopp,  Burke,  and  Cumming  (1972) 
showed  that  the  U.  of  A.  recording  system  with  a  bandpass  of 
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Figure  3.2  Amplitude  and  phase  spectra  of  the  seismic  array 
recording  instrument  used  by  the  University  of  Alberta  in 
1980  . 
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Figure  3.3  Amplitude  and  phase  spectra  of  the  seismometer 
used  by  the  University  of  Alberta  in  1980. 
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0.1  to  50  Hz  coupled  with  a  14-bit  A  to  D  converter  would 
provide  the  wide- f requency  response  and  dynamic  range  for 
high-quality  recording  of  seismic  reflection/  refraction 
signals  from  the  deep  crust,  and  that  the  system  noise  was 
not  a  significant  factor  in  data  analysis. 

Furthermore,  the  U.  of  A.  digital  array  offers 
additional  advantages  by  comparison  with  a  single  channel 
se i smometer : 

1.  An  array  of  detectors  reduces  the  possibility  of 
recording  failure  compared  to  a  system  using  a  single 
se i smometer . 

2.  The  combination  of  individual  records  in  the  array 
improves  the  signal-to-noi se  ratio.  The  s ignal-to-noi se 
ratio  is  increased,  by  a  factor  of  approximately  N0-5,  if 
N  is  the  number  of  outputs  ( B3th , 1 979 ) . 

3.  The  array  of  stations  may  yield  more  information  about 
the  nature  of  the  signals  (e.g  the  apparent  velocity  of 
the  seismic  signals  along  the  profile  line). 


4.  DATA  PROCESSING 


4.1  Comments  on  Raw  Data 

The  refraction  data  for  the  present  study  consist  of 
six  profiles.  Since  we  are  interested  in  only  the  regional 
structure,  no  corrections  were  applied  for  shot  size, 
receiver  elevation,  and  water  depth  for  this  preliminary 
interpretation.  Owing  to  equipment  malfunctions,  5  records 
from  the  U.  of  A.  digital  array  were  found  to  be  unusable, 
and  were  edited  out;  the  amplitudes  of  twelve  records  from 
various  recording  systems  were  also  incorrect  up  to  an  order 
of  magnitude  (Ellis  and  Clowes,  1981,  p.23).  Because  of 

this,  we  normalized  each  trace  to  a  maximum  amplitude  of 
unity  for  plotting  purposes.  Furthermore,  when  the  shots 
were  detonated,  all  receiver  channels  were  not  turned  on  at 
the  same  time.  As  a  result,  the  WWVB  time  is  needed  for  each 
receiver  channel  to  align  them  correctly  with  respect  to  the 
shot  instant. 

In  general,  some  of  the  features  present  in  the  raw 
refraction  data  are: 

1.  The  amplitudes  of  the  nearest  seismic  trace  are 
significantly  larger  than  the  most  distance  ones. 

2.  Beyond  150  km  from  the  shot,  the  amplitudes  of  first 
arrivals  become  small  and  not  distinct  from  the  noise 
background . 

3.  Because  of  high  noise  levels  at  large  distances  and  a 
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small  amplitude  of  Pn ,  only  the  profile  FI  shows  the 
first  arrivals  which  may  be  identified  as  a  Pn  phase. 

4.  A  discontinuity  in  the  travel-time  curve  is  observed  on 
profile  FI  at  a  distance  of  about  210  km  from  the  shot 
point,  but  the  discontinuity  of  the  reversed  profile  N2 
is  not  observed  because  of  the  high  noise  levels. 


4.2  Spectral  Analysis 
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time-domain  measurements  of  the  seismic  trace  into  frequency 
domain  measurements,  thus  the  amplitude  for  each  frequency 
component  can  be  evaluated  through  the  Fourier 
transformation.  A  power  spectrum  is  a  measure  of  the  square 
of  the  amplitude  for  each  frequency.  Knowledge  of  frequency, 
amplitude,  and  phase  characteristics  is  extremely  important 
in  designing  digital  processing  techniques  to  improve 
seismic  data.  Bath  (1974)  provided  a  comprehensive  summary 
of  theory  and  applications  of  spectral  analysis  in  various 
geophysical  desciplines. 

The  finite  length  of  seismic  data  represents  the 
truncation  of  the  infinite  data  by  a  rectangular  function. 
The  discontinuity  between  the  beginning  and  the  end  of  data 
creates  undesirable  side-lobe  effects  in  the  frequency 
domain,  thus  the  power  spectrum  of  the  truncated  data  may  be 


negative  at  some  frequencies.  A  Daniel  window,  which 
represents  the  closest  inverse  of  a  rectangular  window  in 
the  frequency  domain,  is  useful  to  apply  on  data  to 
eliminate  the  undesirable  effect  of  the  truncated  data  and 
to  reduce  the  variance  of  the  power  estimate.  Kanasewich 
(1975,  p.96-105)  presented  a  detailed  discussion  of  the 
theory  and  the  computational  procedure  of  Daniel  power 
estimate  and  he  recommended  the  use  of  the  Daniel  spectral 
estimate  when  the  data  set  had  100  to  4000  samples. 

In  the  time  domain  (or  lag  domain  if  we  calculate  the 
power  spectra  with  the  autocovariance  function)  the  Daniel 
window  is  defined  as  : 

W(k)=sin(ak/m)/ ( ak/m)  k=0,1,2,....n-1 

where : 
a  =  3.1416, 

n  =  the  number  of  observations, 
m  =  the  width  of  the  window. 

The  parameter  m  controls  the  resolution  of  the  spectral 
estimate;  the  greater  the  value  of  m,  the  greater  will  be 
the  resolution.  However,  the  greater  the  value  of  m,  the 
greater  will  be  the  variance  of  the  spectral  estimate  for  a 
given  value  of  n.  Usually  m  is  chosen  to  be  about  10  to  20  % 
of  n . 

A  computer  program  written  by  Robinson  and  Siliva(1979) 
was  used  to  compute  the  Daniel  power  spectrum  in  this  study. 
The  program  first  computes  the  autocorrelation  of  the  data, 
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weights  the  autocorrelation  values  by  the  Daniel  window,  and 
then  collapses  all  values  like  an  accordian  to  length  m  (the 
width  of  the  window).  The  Fourier  transform  converts  the 
collapsed  autocorrelation  values  into  the  frequency  domain, 
thus  the  power  spectrum  is  estimated  at  equi-spaced 
frequencies  between  zero  and  the  Nyquist  frequency. 

Fig. 4.1  to  4.3  show  the  power  spectral  estimates  at 
various  locations  of  increasing  distances  from  the  shot  A1. 
These  figures  indicate  that  the  high  frequency  components  of 
the  seismic  signals  decrease  as  the  distance  increases.  This 
implies  that  the  high  frequency  components  are  filtered  out 
by  the  earth  as  the  seismic  waves  travel  greater  distances. 
In  Fig.  4.1  to  4.3  it  may  be  observed  that  the  upper  limit 
of  significant  power  decreases  from  about  20  Hz  at  a 
distance  of  25  km  to  about  11  Hz  at  152  km.  On  the  whole, 
the  frequency  content  of  'useful'  signal  energy  is 
restricted  to  a  range  of  5  to  20  Hz.  A  bandpass  filter  of 
5-20  Hz  was  applied  to  records  from  profiles  A1,  A2 ,  N1,  and 
F2.  However,  for  the  profile  FI  and  N2,  where  the  distances 
between  the  shots  and  the  recording  stations  were  over  100 
km,  the  records  were  filtered  with  a  bandpass  of  2  to  10  Hz. 

Fig. 4. 4  and  4.5  show  the  dramatic  difference  of 
frequency  contents  in  marine  shots  from  the  profiles  N1  and 
F2.  The  two  shots  were  located  at  the  bottom  of  a  water 
layer  with  depths  of  60  and  110  m  respectively.  The  power 
spectra  exhibit  the  anomaly  of  water-layer  reverberation: 
the  reverberation  of  the  signal  between  the  top  and  bottom 
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Figure  4.1  The  power  spectrum  for  refraction  line  A1.  The 
distance  is  24.6  km  from  the  shot. 
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Figure  4.2  The  power  spectrum  for  refraction  line  A1.  The 
distance  is  41.2  km  from  the  shot. 
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Figure  4.3  The  power  spectrum  for  refraction  line 
distance  is  152.0  km  from  the  shot. 
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Figure  4.4  The  power  spectrum  of  marine  shooting  for  line 
N1.  Note  the  large  amplitude  at  6.0  Hz  due  to  water 
reverberat ion . 
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Figure  4.5  The  power  spectrum  of  marine  shooting  for  line 
F2.  Note  the  large  amplitude  at  3.7  Hz  due  to  water 
reverberation . 
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of  the  water  layer  gives  rise  to  a  sharply  defined  maximum 
in  the  frequency  response.  The  two  peak  frequencies  are  at 
3.7  and  6.0  Hz  respectively.  The  application  of 
deconvolution  filters  is  often  suggested  as  a  mean  of 
removing  reverberation  effects  which  obscure  the  recorded 
signals.  These  filters  might  be  useful  for  some  of  the 
present  data,  particularly  when  the  later  arrivals  are  to  be 
studied . 

On  the  other  hand,  the  power  spectrum  of  the  water 
reverberation  is  useful  to  estimate  the  depth  of  the  water 
layer.  The  relation  between  the  peak  frequency  and  the  water 
depth  is  expressed  (Dobrin,  1976)  by 

f =v/4d 

where : 

f  =  peak  frequency, 
v  =  velocity  of  sound  in  water, 
d  =  water  depth. 

Two  examples  are  used  to  illustrate  that  the  water 
depth  can  be  estimated  by  the  power  spectrum  of  water 
reverberation  using  the  above  formula.  Table  1  summarizes 
the  results  of  the  measured  water  depths  (d,)  (Ellis  and 
Clowes,  1981,  p . 1 2 )  and  water  depths  (d2)  determined  from 
the  power  spectra. 

Comparison  of  the  measured  water  depths  and  depths  by 
power  spectra  from  the  two  seismic  records  (Table  1)  shows 
only  minor  differences.  This  indicates  that  the  depth  of  the 
water  layer  can  be  estimated  from  the  power  spectrum  of  data 
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geophone 

v (m/s ) 

f  (Hz  ) 

d,  (m) 

d  2 

6002  1 

1490 

6.0 

60 

62 

45004  1 

1490 

3.7 

1  10 

100 

Table  1....  Comparison  of  the  results  of  the  measured  water 
depths  to  the  depths  determined  from  power  spectra. 


which  contain  the  signals  of  water  reverberation. 


4.3  Bandpass  Filtering 

A  bandpass  filter  is  a  frequency  filter  designed  to 
pass  signal  frequencies  in  a  particular  band  and  to 
attenuate  all  other  frequencies.  The  purpose  of  applying  a 
bandpass  filter  to  seismic  signals  is  to  extract  useful 
signals  from  noise  background.  The  operation  of  filtering  is 
especially  effective  in  cases  where  signal  and  noise  spectra 
do  not  overlap  over  a  wide  frequency  band.  Two  types  of 
bandpass  filters  are  considered  in  this  study:  eight-pole, 
zero-phase  shift,  recursive  Butterworth  and  Bessel  bandpass 
filters.  The  effectiveness  of  both  filters  is  compared  on 
the  data  containing  the  first  arrival  of  energy,  and  noisy 
data . 

Kanasewich  (1975,  p.175-203)  provided  a  detailed 
discussion  of  the  8-pole  Butterworth  filter.  By  specifying 
the  attenuation  rate  at  a  high  and  low  cutoff  frequencies, 
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the  filter  coefficients  are  evaluated  from  poles  and  zeros 
of  the  Butterworth  polynomial  in  a  Z-plane  representation. 
The  input  data  are  then  convolved  with  filter  coefficients 
recursively  in  both  forward  and  backward  directions  to 
achieve  zero  phase  shift  filtering. 

The  design  of  a  8-pole  Bessel  filter  is  based  on  the 
Bessel  polynominal .  The  derivation  of  filter  coefficients  is 
carried  out  in  a  similar  manner  to  those  for  the  Butterworth 
filter.  A  good  discussion  of  the  Bessel  filter  is  given  by 
Lubk in  (1970). 

Comparison  of  the  Bessel  and  Butterworth  filters  shows 
the  following  differences: 

1.  The  Butterworth  filter  has  a  much  faster  transition 
region  from  passband  to  stopband,  and  has  a  much  better 
amplitude  characteristic  than  that  of  the  Bessel  filter. 

2.  In  the  time  domain,  the  Bessel  filter  has  a  only  a  tenth 
the  overshoot  of  the  Butterworth  for  filters  up  to  five 
poles,  and  gets  better  as  more  poles  are  added,  while 
the  Butterworth  filter  gets  worse. 

Both  filters  were  applied  to  the  data  that  contained 
the  first  arrival  of  energy.  Fig.  4.6  to  4.8  show  the 
application  of  both  filters  on  the  first  arrival  over 
various  frequency  bands.  The  Bessel  filter  has  little 
distortion  on  the  first  arrival  in  a  range  of  frequency 
bands,  but  the  Butterworth  filter  introduces  a  significant 
offset  on  the  first  arrival,  as  the  frequency  band  gets 
narrower.  For  example,  at  the  frequency  band  5-12  Hz,  the 
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8  POLES  BUTTERWGRTH  BANDPASS  FILTER 


8  POLES  BESSEL  BANDPASS  FILTER  1-20  HZ 


Figure  4.6  Comparison  of  Butterworth  and  Bessel  bandpass 
filters  on  the  first  arrival  in  the  frequency  band  of  1-20 
Hz. 


32 


TIME  (SEC) 


8  POLES  BUTTERWGRTH  BANDPASS  FILTER 


8  POLES  BESSEL  BANDPASS  FILTER  5-20  HZ 


Figure  4.7  Comparison  of  Butterworth  and  Bessel  bandpass 
filters  on  the  first  arrival  in  the  frequency  band  of  5*20 
Hz . 
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TIME  (SEC) 


8  POLES  BUTTERNORTH  BRNDPRSS  FILTER 


8  POLES  BESSEL  BRNDPRSS  FILTER  5-12  HZ 


Figure  4.8  Comparison  of  Butterworth  and  Bessel  bandpass 
filters  on  the  first  arrival  in  the  frequency  band  of  5-12 
Hz. 
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Butterworth  filter  offsets  the  first  arrival  time  by  0.2 
sec.  This  indicates  that  the  Bessel  filter  is  more  suitable 
to  apply  on  the  data  that  contains  first  arrival  energy,  for 
the  Bessel  has  a  rather  smooth  transition  region  from  the 
passband  to  stopband  and  tapers  off  the  input  data  more 
gradually,  thus  avoiding  distortion  of  the  signals. 

Both  filters  were  then  tested  on  a  noisy  record.  Fig. 
4.9  and  4.10  display  the  raw  data  and  the  data  filtered  by 
both  filters.  Comparison  of  the  two  filtered  records  shows 
the  following  features: 

1.  The  Butterworth  filter  clearly  suppresses  the 
undesirable  high  and  low  frequency  noise  outside  the 
passband . 

2.  The  Bessel  filter  is  not  able  to  suppress  completely 
both  high  and  low  frequency  noise  outside  the  passband. 
For  example,  the  low  frequency  noise  about  1.3  Hz  (trace 
3)  is  still  present  in  the  record. 

3.  Both  filtered  records  have  improved  the  correlation  of 
seismic  traces  considerably. 

4.  The  appearances  of  both  filtered  records  are  very 
similar  to  each  other. 

The  Butterworth  filter,  in  general,  has  better 
performance  than  that  of  the  Bessel  filter  in  extracting 
signals  from  noise  background.  The  rapid  transition  from 
passband  to  stopband  of  the  Butterworth  is  quite  good  to 
obtain  a  specific  band  of  frequencies,  and  efficient  to 
attenuate  other  frequencies  that  are  outside  the  passband. 
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TIME  (SEC) 


REFRACTION  DATA  L  I NE--FUCA  (A2) 
BANDPASS  FILTER  BESSEL  2  TO  10  HZ 


TIME  (SEC) 


Figure  4.9  The  raw  refraction  record  and  the  corresponding 
record  filtered  by  a  Bessel  bandpass  filter. 
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TIME  (SEC) 


REFRACTION  DATA  L  I  NE--FUCA  ( A2 ) 
BANDPASS  FILTER  BUTTERHGRTH  2  TO  10  HZ 


TIME  (SEC) 


Figure  4.10  The  raw  refraction  record  and  the  corresponding 
record  filtered  by  a  Butterworth  bandpass  filter. 
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However,  since  the  travel  times  of  first  arrivals  are  our 
main  interest,  the  minimum  distortion  on  the  first  arrivals 
introduced  by  the  Bessel  filter  is  more  desirable.  Thus  we 
decided  that  the  Bessel  filter  would  be  used  on  all  seismic 
records . 


It  should  be  noted  that  in  the  form  in  which  the  Bessel 
filter  was  used  in  this  work,  the  effective  bandpass  is 
considerably  wider  than  that  of  the  Butterworth  filter  with 
the  same  specified  frequencies.  Thus  the  poorer  low 
frequency  rejection  of  the  Bessel  filter  is  a  function  of 
larger  bandwidth  and  could  be  improved  by  narrowing  the 
nominal  bandwidth. 


5.  VELOCITY-DEPTH  DETERMINATION  FROM  REFRACTION  DATA 


5.1  Comments  on  Refracted  Waves  from  the  Crust 

Since  attenuation  and  absorption  of  elastic  waves 
through  the  earth  cause  the  amplitudes  of  elastic  waves  to 
decrease  with  increasing  distance,  study  of  the  amplitude  is 
essential  to  understand  how  the  head  wave  propagates  through 
the  earth.  Grant  and  West  (1965,  p . 181 )  summarized  three 
important  facts  for  the  propagation  of  head  waves  through 
the  earth:  1.  the  head-wave  amplitude  diminishes  inversely 
as  the  square  of  the  distance  from  the  source;  2.  since  the 
absorption  coefficient  is  directly  proportional  to 
frequency,  the  head-wave  amplitude  also  diminishes  with 
frequency  with  the  increase  of  travel  time;  3.  the  energy  in 
the  head  wave  becomes  maximum  at  the  critical  distance  (the 
point  where  the  refracted  wave  meets  the  reflected 
wavefront ) . 

Based  on  travel-time  analysis,  three  dominant  seismic 
waves  are  observed  from  the  surface  down  to  the  upper 
mantle : 

1.  Pg  is  usually  assigned  to  be  a  refracted  wave  in  the 
upper  crust,  which  has  a  typical  velocity  of  5.5  to  6.1 
km/s.  The  Pg  phase  is  usually  observed  up  to  a  distance 
of  60-100  km.  In  general,  the  amplitude  of  the  Pg  phase 
decays  inversely  as  the  square  of  the  distance  from  the 
source;  beyond  100  km  its  amplitude  becomes  very  small, 
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and  sometimes  is  not  even  observable. 

2.  The  next  prominent  phase,  which  travels  in  the  lower 
crust,  is  designated  as  P+ .  However,  this  prominent 
phase  is  only  observed  in  some  regions  (Berry,  1973). 
The  absence  of  P+  phase  is  probably  due  to  lateral 
heterogeneities  in  lower  crustal  layer.  The  typical 
velocity  is  between  6.8  to  7.1  km/s. 

3.  Pn  is  generally  considered  to  be  a  head  wave  which 

travels  directly  beneath  the  Moho  discontinuity.  In 
general,  the  Pn  phase  is  observed  to  be  a  first  arrival 
at  distances  beyond  120-200  km.  Because  of  the 

attenuation  and  absorption  of  elastic  waves  through  the 
earth  the  amplitude  of  Pn  decays  rapidly  beyond  the 
critical  distance.  The  typical  Pn  phase  velocity  in 
British  Columbia  is  between  7.8  to  8.1  km/s. 


5.2  Record  Representation 

All  seismic  sections  were  represented  on  a  reduced 
travel-time  diagram,  where 


T  =  t-D/V 

wi  th : 

T  =  reduced  travel  time, 
t  =  travel  time  of  seismic  arrivals, 
D  =  source-receiver  distance, 

V  =  reduction  velocity. 
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This  kind  of  display  of  data  gives  much  better 
resolution  on  a  given  graph  when  large  recording  distances 
are  involved,  since,  when  the  reduction  velocity  is  close  to 
the  phase  velocity  of  particular  events  of  interest,  these 
are  more  or  less  horizontal  on  the  diagram.  This 
reduced-time  representation  offers  the  advantage  of 
representing  more  time  data  within  a  given  graph,  or 
conversely  allowing  an  expanded  time  scale  to  be  used.  Since 
the  Moho  discontinuity  was  of  a  particular  interest  in  this 
study,  a  reduction  velocity  of  8.0  km/s  was  applied  on  all 
record  sections. 


5.3  Correlation  of  Refraction  Arrivals 

Correlation  of  first  arrivals  from  trace  to  trace  is 
essential  to  determine  the  apparent  phase  velocity  of  a 
particular  travel-time  branch  in  the  refraction 
interpretation.  In  general,  the  attenuation  and  absorption 
of  elastic  waves  in  rock  cause  a  progressive  lowering  of 
apparent  frequency  of  seismic  events  with  increasing 
distance  of  travel  through  the  earth,  thus  the  shape  of  the 
waveform  changes  as  a  function  of  distance.  In  addition, 
when  noise  is  present  in  some  records,  it  is  difficult  to 
observe  the  exact  instant  of  onset  of  the  first  break 
energy.  Thus,  it  is  common  practice  not  to  use  the  beginning 
of  the  event  for  correlation  purposes,  but  the  pronounced 
maximum  or  minimum  amplitude  following  the  first  arrivals  of 


4  1 


energy.  The  correlation  of  convenient  peaks  or  troughs  of 
the  seismic  signals  which  have  similar  appearances  provides 
a  more  reliable  estimate  of  apparent  velocity  of  the 
travel-time  branch.  In  order  to  determine  the  actual  arrival 
times,  we  need  to  make  a  phase  correction  from  the  picked 
value  of  the  peak  or  trough,  back  to  the  estimated  beginning 
of  the  refraction  pulse. 


5.4  Plane-Layered  Interpretation 

First  arrivals  were  picked  from  all  record  sections, 
except  the  last  two-thirds  of  the  first  arrivals  on  profile 
N2.  Linear  segments  of  travel  times  were  fitted  in  a 
least-square  program  to  determine  velocities  and  intercept 
times  of  the  refraction  waves  that  travel  at  different 
depths.  The  determination  of  the  near  surface  velocity  was 
not  well  defined  due  to  the  relatively  large  separation  of 
recording  sites.  The  travel  time  to  the  first  recording  site 
from  four  shot  locations  yields  a  range  of  near  surface 
velocities  from  14.1  km/s  to  4.57  km/s  (Table  2).  The  one 
extremely  high  apparent  near  surface  velocity  (14.1  km/s  ) 
is  definitely  due  to  an  error  of  the  start  time,  therefore, 
it  is  reasonable  to  ignore  this  high  surface  velocity.  The 
velocity  from  the  shot  A1  and  F2  are  close  to  an  average 
velocity  of  5.3  km/s.  This  average  velocity  is  considered  to 
be  an  upper  limit  of  near  surface  velocity  and  can  be  used 
to  compute  the  necessary  models  for  the  preliminary 
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refraction 

line 

geophone 

station 

di stance 
km 

t  ime 
sec 

surface 
veloc i ty 
km/s 

A 1 

4500  1  1 

3.308 

0.626 

5.28 

A2 

10051 

3.671 

0.803 

4.57 

F2 

45004 1 

14 . 743 

2.766 

5.33 

N  1 

1  002  1 

2.036 

0.144 

14.11 

Table  2....  Near  surface  velocities  from  seismic  profiles  of 
A  1  ,  A2 ,  F2 ,  and  N 1  . 


refraction  interpretation. 

The  depth  to  the  refractor  was  calculated  from  the 
intercept  time  using  a  horizontal  plane-layered  model.  A 
computer  program  listed  in  the  Appendix  allows  the 
calculation  of  the  models  for  n-horizontal  layers. 

Central  Profiles 

Fig. 5.1  shows  the  Argo  refraction  profile  ( A 1 ) .  This 
profile  consists  of  a  near  surface  layer  and  two  distinct 
deeper  layers.  The  velocities  of  these  layers  are  6.35  and 
7.02  km/s.  The  corresponding  thicknesses  are  0.51  km  for  the 
surface  layer  and  10.5  km  for  the  upper  crust.  The  noise 
levels  in  this  section  are  not  significant;  the  first 
arrivals  have  relatively  large  amplitudes  and  correlate  very 
well  throughout  the  whole  section. 


REFRACTION  LINE  ARGO  ( A 1  ) 
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Figure  5.1  Vertical  component  reduced  section  of  line  A1. 
Small  ticks  on  traces  indicate  first  arrival  picks. 
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Fig. 5. 2  displays  the  record  section  of  profile  A2 .  This 
section  is  similar  to  section  A1.  The  first  arrivals  from 
both  lower  and  upper  crust  are  clear  and  well  defined.  Their 
velocities  are  6.39  and  6.99  km/s,  corresponding  to 
thicknesses  of  0.58  and  13.1  km  respectively.  Comparison  of 
crustal  velocities  with  profile  A1  shows  that  both 
velocities  are  in  excellent  agreement,  but  the  thickness  of 
the  upper  crust  has  a  rather  large  discrepancy  between  two 
profiles.  This  discrepancy  is  up  to  2.6  km.  This  implies 
that  the  crustal  structure  beneath  the  shot  A  is  not  a 
simple  horizontal  model,  and  that  lateral  variations  of 
structure  may  be  present. 

Nero  Profiles 

Fig. 5. 3  and  5.4  are  profiles  of  the  observations  from 
the  northern  end  along  the  structural  axis  of  Vancouver 
Island  to  the  southern  end  of  the  island.  Fig.  5.3  is  the 
reversed  profile  of  Fig.  5.1.  The  first  arrivals  are  clear 
and  coherent  up  to  a  distance  of  210  km  from  the  shot.  The 


data  were  well 

f i t ted 

into  two 

linear 

segments 

which 

correspond  to 

the 

veloc i t ies , 

6.57 

and  7.04 

km/s 

respectively.  The 

very 

high  noise 

levels  and 

small 

amplitudes  of  refracted  first  arrivals  made  the  correlation 
beyond  210  km  very  difficult.  A  narrow  bandpass  (2  Hz  to  5 
Hz)  was  applied  to  the  data  in  attempt  to  recover  signals 
from  the  noisy  background,  but  the  low  s i gnal-to-noi se  ratio 
made  travel-time  picks  ambiguous  so  that  the  determination 


REFRRC T I  ON  LINE  RRGO  ( R2  ) 
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Figure  5.2  Vertical  component  reduced  section  of  line  A2. 
Small  ticks  on  traces  indicate  first  arrival  picks. 
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Figure.  5.3  Vertical . component  reduced  section  of  line  N1 
Small  ticks  on  traces  indicate  first  arrival  picks. 
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Figure  5.4  Vertical  component  reduced  section  of  line 
N2. Small  ticks  on  traces  indicate  first  arrival  picks. 
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of  Pn  phase  velocity  was  not  possible. 

Fuca  Profiles 

Fig. 5. 5  and  5.6  represent  the  reversed  section  of  the 
Nero  profile.  The  profile  F2  is  similar  to  N1.  The  first 

arrivals  of  6.52  km/s  continues  up  to  100  km,  but  the 

arrivals  from  the  lower  crust  are  scattered.  They  were 

fitted  approximately  with  a  6.95  km/s  line  segment.  A 

discontinuity  of  travel-time  curve  is  observed  about  210  km 
from  the  shot  on  the  profile  FI.  This  discontinuity  may 
indicate  the  existence  of  a  low  velocity  zone.  If  such  a  low 
velocity  zone  exists  in  this  region,  it  is  difficult  to 
establish  the  depth  to  the  Moho  from  the  refraction  data. 
Arrivals  identified  as  representing  a  Pn  phase  are  observed 
on  the  profile  FI  between  the  distance  220  to  290  km,  but 
their  amplitudes  are  quite  small.  In  order  to  enhance  the 
certainty  of  the  Pn  arrivals,  we  enlarged  the  amplitude  of 
Pn  arrivals  and  reduced  the  background  noise  considerably  by 
applying  a  Butterworth  filter  with  a  bandpass  of  1  to  2.5  Hz 
on  the  profile  of  FI.  The  Pn  arrivals  are  clear  and  coherent 
in  Fig. 5. 7.  The  apparent  velocity  was  found  to  be  7.95  km/s. 
By  assuming  plane  horizontal  layers,  the  depth  for  the  Moho 
determined  from  intercept  times  is  35  km  (  Fig.  5.6).  Later 
arrivals  in  Fig.  5.7,  which  may  be  reflections  from  the 
Moho,  appear  fairly  clear  and  coherent,  with  rather  large 
amplitudes.  Another  feature  noteworthy  in  Fig.  5.7  is  the 
line  C-C,  postulated  to  be  the  secondary  arrivals  from  a 
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Figure  5.5  Vertical  component  reduced  section  of  line 
Small  ticks  on  traces  indicate  first  arrival  picks. 
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Figure  5.6  Vertical  component  reduced  section  of  line  FI. 
Small  ticks  on  traces  indicate  first  arrival  picks. 
PmP=ref lec t ions  from  the  Moho. 
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Figure  5.7  Vertical  component  reduced  section  of  line  F2. 
Small  ticks  on  traces  indicate  first  arrival  picks. 
PmP=ref lect ions  from  the  Moho;  CC=post ula ted  secondary  Pn 
arrivals  from  a  sub-Moho  refractor. 
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sub-Moho  refractor.  These  secondary  arrivals  have  fairly 
large  amplitudes  and  an  apparent  velocity  of  about  8.3  km/s. 
Because  of  the  truncation  of  the  profile  at  330  km,  first 
break  arrivals  are  not  recorded  for  this  later  branch  of  Pn 
arrivals.  Thus  the  interpretation  of  this  later  branch  is 
not  attempted. 


5.5  Wiechert-Herglotz-Bateman  (WHB)  Inversion  of  the 
Travel-Time  Curves 

Although  the  simplified  approach  of  a  horizontal  layer 
model  provides  useful  information  of  an  approximate 
velocity-depth  function,  it  requires  the  calculation  of 
thicknesses  by  a  series  of  possible  layers  with  a  constant 
velocity  for  each  layer.  In  reality,  velocity  usually 
increases  with  depth.  Another  approach,  which  makes  use  of 
velocity  increasing  continuously  with  depth,  may  represent  a 
more  realistic  model  of  the  earth.  This  approach  is  known  as 
inversion  of  travel-time  curve  by  WHB  integral.  Grant  and 
West  (1965,  p. 138—141 )  provided  a  detailed  analysis  of  this 
approach . 

We  assume  that  the  elastic  waves  travel  through  a 
succession  of  horizontal  layers  of  different  elastic  media, 
and  the  velocity  increases  as  a  function  of  depth  only.  The 
WHB  integral  for  plane  horizontal  layers  may  be  shown  to  be: 
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where  : 

D  =  the  distance  of  a  ray  path  that  returns  to  the  surface 
of  the  ground, 

V  =  the  velocity  of  the  ray  at  the  maximum  depth  of 
penetration , 

dt/dx  =  the  derivative  of  the  time-distance  curve. 

Before  the  application  of  the  integral,  the  travel-time 
data,  excluding  the  data  for  the  Pn  arrivals,  were  reduced 
to  the  base  of  the  surface  layer  by  subtracting  the  times 
that  were  calculated  from  the  intercept  times  of  the 
plane-layered  model.  This  removes  the  effect  of  the  near 
surface  layer.  The  travel-time  data  were  then  fitted  by  a 
polynomial  of  arbitrary  degree.  The  use  of  the  polynomial 
has  certain  advantages  in  analysing  the  data:  1.  we  can 
establish  a  concise  mathematical  formula  to  express  the 
relation  between  two  variables  (in  our  case,  time  and 
distance);  2.  although  the  degree  of  the  polynomial  to  be 
fitted  to  the  data  may  not  be  known  in  advance,  we  change 
the  order  of  the  polynomial  until  we  reach  an  order  for 
which  we  get  a  satisfactory  result  of  fitting  the  data;  3. 
the  polynomial  is  a  continuous  function  and  is 
differentiable  anywhere  within  its  domain,  it  is  useful,  for 
instance,  for  interpolation  purposes  or  for  smoothing  the 
original  data  prior  to  further  analysis. 

In  general,  the  degree  of  polynomial  is  not  restricted. 
However,  one  condition  is  required  in  order  to  evaluate  the 
WHB  integral,  that  is:  V  >  dx/dt .  We  found  that  the  data  of 
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sections  A1,  A2,  and  Hero  fitted  with  a  polynomial  of  second 
degree;  the  data  of  section  Fuca,  with  the  polynomial  up  to 
fourth  degree.  The  coefficients  of  the  polynomial,  in  turn, 
permit  the  evaluation  of  the  WHB  integral  by  determining  the 
derivative  of  the  time-distance  curve  at  any  arbitrary 
distance.  A  computer  program  listed  in  the  Appendix  was 
written  for  the  inversion  of  the  travel-time  data  by  this 
method . 

Fig. 5. 8  to  5.11  illustrate  the  application  of  the  WHB 
integral  to  the  time-distance  curves  of  sections  A1,  A2, 
Nero,  and  Fuca.  The  velocity-depth  functions  of  the  four 
seismic  sections  are  quite  similar  to  a  depth  of  6  km;  the 
velocity  beneath  the  surface  layer  of  5.3  km/s  gradually 
increases  with  depth  from  the  initial  velocity  of  about  6.2 
km/s  to  a  velocity  6.5  km/s  at  a  depth  of  6  km.  Below  this 
depth  the  velocity-depth  curve  of  profiles  A2  and  Fuca  have 
similar  slopes,  but  their  slopes  are  steeper  than  those  from 
the  profiles  Nero  and  A1.  As  a  result,  the  depth  to  any 
given  velocity  refractor  from  the  profiles  A2  and  Fuca  is 
greater  than  that  from  profiles  Nero  and  A1.  On  the  other 
hand,  the  velocity-depth  curve  of  profile  A1  has  a  minimum 
slope  which,  in  turn,  yields  a  minimum  depth  to  a  given 
velocity  refractor.  For  example,  the  depth  to  a  velocity  of 
7.0  km/s  fluctuates  from  a  maximum  of  21.4  km  from  the 
profile  A2  to  a  minimum  of  16.5  km  from  the  profile  A1,  with 
a  intermediate  value  of  17.3  km  from  the  profile  Nero.  This 
implies  that  the  ray  path  travelled  from  A  to  F  is  quite 
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REFRACTION  LINE  ARGO  (All 


DEPTH  KM 

Figure  5.8  Time-distance  curve:  cross=  raw  data;  solid  line= 
data  fitted  by  a  polynomial,  and  the  velocity-depth  curve  by 
the  WHB  inversion  of  time-distance  curve  beneath  the  surface 
layer . 
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REFRACTION  LINE  ARCO  (A2) 


DEPTH  KM 

Figure  5.9  Time-distance  curve:  cross=  raw  data;  solid  line= 
data  fitted  by  a  polynomial,  and  the  velocity-depth  curve  by 
the  WHB  inversion  of  time-distance  curve  beneath  the  surface 
layer . 
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REFRACTION  LINE  NERO  (N1-N2) 


DEPTH  KN 


Figure  5.10  Time-distance  curve:  cross=  raw  data;  solid 
line=  data  fitted  by  a  polynomial,  and  the  velocity-depth 
curve  by  the  WHB  inversion  of  time-distance  curve  beneath 
the  surface  layer. 
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REFRACTION  LINE  FUCA  (F1-F2) 


DEPTH  KM 


Figure  5.11  Time-distance  curve:  cross=  raw  data;  solid 
line=  data  fitted  by  a  polynomial,  and  the  velocity-depth 
curve  by  the  WHB  inversion  of  time-distance  curve  beneath 
the  surface  layer. 


■ 
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different  from  the  ray  path  that  travels  from  A  to  N.  One 
possible  explanation  of  the  discrepancy  of  travel  times  in 
the  southerly  and  northerly  directions  is  the  existence  of  a 
zone  of  anomalous  structure  south  of  the  shot  point  A.  If 
this  hypothesis  is  true,  it  might  account  for  some  of  the 
discrepancies  in  the  travel-time  curves.  The  extent  of  the 
anomalous  zone  cannot  be  estimated  precisely  from  the 
refraction  data  of  the  present  study. 


5.6  Synthetic  Seismograms 

The  interpretation  based  on  travel-time  analysis  alone 
is  not  always  unique,  and  a  comparison  of  theoretical  and 
experimental  travel-time  curves  can  be  useful  to  avoid 
mistakes  in  the  solution  of  an  inverse  problem.  The 
theoretical  travel  times  can  be  approximated  by  modified 
asymptotic  wave  methods  (Cervery  and  Zahradni k , 1 972 ; 
Choi , 1978) . 

The  method  used  to  compute  time-distance  curves  in  the 
present  study  is  based  on  a  modified  Airy  function  solution 
derived  by  Choi(1978).  His  approach  was  to  derive  a  integral 
solution  for  an  arbitrary  ray  in  a  layered,  vertically 
inhomogeneous  elastic  medium,  using  a  modified  version  of  a 
third-order  saddle  point  approximation.  This  result  was  then 
reduced  to  the  form  of  a  modified  Airy  function  solution  for 
the  computation  of  a  synthetic  seismogram.  For  a  detailed 
analysis  the  reader  is  referred  to  Choi(1978). 
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Since  we  anticipate  lateral  velocity  variations  and  the 
existence  of  a  low  velocity  zone  as  suggested  by  the  first 
arrival  data,  we  do  not  expect  that  all  the  features  on  the 
observed  record  section  can  be  matched  by  a  modelling 
procedure  using  synthetic  seismograms.  In  general,  the 
theoretical  computations  provide  an  approximate 
velocity-depth  structure  which  is  at  least  consistent  with 
the  observed  travel-time  and  amplitude  data. 

The  construction  of  the  synthetic  seismogram  was 
initially  based  on  the  velocity-depth  function  of  the  Fuca 
profiles  (F1-F2)  which  was  derived  from  the  WHB  inversion  of 
travel  times.  The  parameters,  velocity  as  well  as  depth,  of 
the  synthetic  seismogram  were  adjusted  to  fit  with  the 
observed  data  of  the  Fuca  profiles.  The  data  to  be  fitted  by 
the  synthetic  seismogram  include  the  refractions  from  the 
crust  and  from  the  mantle,  and  reflections  from  the  mantle. 
In  order  to  improve  the  fit  of  the  observed  data,  a  low 
velocity  zone  was  required  to  keep  the  reflected  travel 
times  correct  at  larger  distances  and  to  generate  the  large 
amplitudes  of  the  possible  Moho  reflections  described  above. 
The  Moho  reflections  from  the  synthetic  seismogram  fitted 
very  well  to  observed  data  of  the  Fuca  profile  (in  Fig. 
5.6).  The  critical  distance  from  the  synthetic  calculation 
was  at  140  km.  Unfortunately,  the  synthetic  seismogram 
program  was  not  designed  to  compute  the  theoretical  travel 
times  of  the  Pn  phase  branch,  however,  the  theoretical  Pn 
travel  times  may  be  estimated  from  the  critical  distance  as 
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determined  from  the  synthetic  seismogram  with  the 
assumptions  of  a  horizontal  Moho  discontinuity  and  a  mantle 
velocity  of  8.0  km/s  (rounded  off  from  the  observed  velocity 
of  7.95  km/s).  The  observed  Pn  arrival  times  are 
approximately  one  second  earlier  than  those  predicted  by  the 
above  model.  We  then  varied  the  velocity  and  the  thickness 
of  the  low  velocity  zone  in  an  attempt  to  improve  the  fit  of 
the  Pn  arrivals.  This  attempt  failed.  We  could  fit  either 
the  reflected  or  refracted  branch  from  the  Moho 
discontinuity  to  the  observed  data,  but  were  unsuccessful  in 
fitting  both  simultaneously. 

Since  we  anticipate  lateral  velocity  variations  beneath 
Vancouver  Island,  the  assumptions  of  a  horizontal  Moho 
discontinuity  and  a  constant  velocity  over  100  km  may  not  be 
valid  for  the  computation  of  the  theoretical  travel  times  of 
Pn  arrivals.  Significant  lateral  velocity  variations  may 
account  for  the  early  Pn  arrivals  of  the  observed  data. 
Another  plausible  reason  is  dip  on  the  Moho  discontinuity.  A 
few  degrees  of  dip  over  100  km  is  sufficiently  large  to 
produce  time  shifts  of  the  magnitude  required  to  bring 
reflected  and  refracted  phases  into  an  agreement.  Since  the 
Pn  arrivals  in  the  Nero  profile  (reversed  profile)  can  not 
be  identified,  the  determination  of  the  dip  angle  is  not 
possible.  If  the  Moho  discontinuity  has  dip,  the  observation 
of  the  earlier  Pn  arrivals  requires  the  waves  to  propagate 
in  the  updip  direction,  i.e.  the  Moho  interface  is  dipping 
from  the  northern  end  towards  the  southern  end  of  the 
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island.  We  may  estimate  the  magnitude  of  time  shifts  caused 
by  the  effect  of  dip  angle  by  considering  a  dipping 
two-layer  structure.  We  assume  the  upper  layer  velocity  to 
be  6.5  km/s ,  and  lower  layer  velocity  to  be  8.0  km/s.  Only 
about  3.0  degrees  of  dip  is  required  to  produce 
approximately  1.0  sec  to  1.5  sec  time-shift  difference 
between  the  two-layer  horizontal  and  dipping  model,  at 
distances  200  km  to  300  km  from  the  shot.  If  the  Moho 
beneath  Vancouver  Island  has  3  degrees  of  dip,  its  true 
velocity  will  become  approximately  7.7  km/s. 

The  model  that  produced  the  best  fit  to  the  observed 
refractions  from  the  crust  and  reflections  from  the  mantle 
for  the  Fuca  profiles  has  following  features:  1.  the 
velocity  increases  continuously  to  a  depth  of  27  km;  2. 
below  this  depth  the  velocity  gradually  decreases  to  6.2 
km/s  at  the  depth  of  38  km;  3.  the  low  velocity  zone  stays 
constant  at  6.2  km/s  for  a  thickness  of  6  km  and  sharply 


jumps  to  a  velocity  of  8.0 

km/ s .  Fig  .  5 .  1  2 

shows  the 

ray 

diagram  and  travel-time 

curves  of 

this 

model . 

The 

travel-time  curves  generated 

by  this  model 

agree 

very 

well 

with  the  observed  data.  This  model  provides  an  alternative 
interpretation  which  may  be  compared  to  the  plane-layered 
model  derived  earlier. 
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Figure  5.12  Synthetic  ray  diagram  and  corresponding 
travel-time  curves  for  the  Fuca  profile:  sol i d=obse rved , 
tr iangle=computed;  branch  1=  refractions  from  the  crust, 
branch  2=  reflections  from  the  mantle. 
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5.7  Particle  Motion  Analysis 

Particle  motion  diagrams  may  be  constructed  by  plotting 
the  radial  and  vertical  seismometer  signals  as  a  function  of 


time  for  the  first  arrival  data.  The 

or  ientat ion 

of 

the 

principal  axis 

of 

the  particle 

trajectory 

in 

the 

radial-vertical 

plane 

indicates  the 

apparent  angle 

of 

emergence  of  the  refracted  wave.  Montalbetti  (1969)  and 
Kazmierczak  (1980)  provided  good  discussions  of  particle 
motion  analyses  of  both  compress i onal  and  shear  type 
motions . 

Three  components  of  ground  motion  (vertical, 
north-south  horizontal,  and  east-west  horizontal)  for  all 
shots,  were  recorded  at  a  few  stations  in  Vancouver  Island. 
All  component  data  were  recorded  on  the  Sprengether  and  U. 
of  A.  digital  array  systems.  Unfortunately,  the  information 
from  the  U.  of  A.  digital  array  were  not  applicable  for 
particle  motion  analysis  either  because  of  noisy  data  or 
dead  traces.  Thus,  the  analysis  was  based  only  on  the  data 
recorded  by  the  Sprengether  system  from  four  different  shot 
points . 

The  radial  component  is  obtained  by  rotating  the  two 
horizontal  components  into  a  direction  corresponding  to  a 
great  circle  azimuth  from  source  to  receiver.  The  azimuth  is 
the  angle  measured  clockwise  from  the  north.  Two  methods 
were  used  to  estimated  the  azimuth: 

1.  Theoretical  calculation  of  azimuth  from  the  geometry  of 
both  source  and  receiver  coordinates. 
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2.  Direct  measurement  of  azimuth  from  the  plot  of  the 
north-south  and  east-west  horizontal  components. 

The  maximum  difference  of  the  azimuth  obtained  from  both 
methods  is  up  to  20  degrees.  We  decided  to  use  the  azimuth 
measured  from  the  plot  of  north-south  and  east-west 
components  to  transform  to  radial  motion  because  we  believed 
that  the  orientation  of  the  detectors  in  the  field  was  only 
approximately  in  the  north-south  direction,  and  that  the 
measured  azimuth  was  more  reliable  than  the  azimuth  computed 
by  the  theoretical  calculation. 

Fig. 5. 12  to  5.16  show  the  first  arrival  events  of  three 
horizontal  components  and  one  vertical  component  for 
different  shot  locations.  It  is  obvious  that  P-wave  type 
motion  is  present  in  all  records.  The  apparent  angle  of 
emergence  measured  from  the  horizontal  varies  from  31  to  51 
degrees . 

Since  the  tilt  angle  of  the  principal  axis  of  the 
particle  trajectory  is  equivalent  to  the  apparent  angle  of 
emergence  of  the  refracted  wave,  we  could  make  use  of 
amplitudes  Av  and  Ah  of  the  vertical  and  horizontal  ground 
displacement  to  measure  the  apparent  angle  of  emergence 
which  is  defined  as 

a  =  tan " 1  ( Av/Ah ) 

where:  a  is  the  apparent  angle  of  emergence  measured  with 
horizontal  boundary.  Walker  (1919)  derived  a  relation 
between  the  apparent  and  true  angle  of  emergence  in  the  form 


2cos2/?  =  (  1 -y  ) /(  0 . 5-y  )  (  1 -Sina  ) 


Mcmit. 
VERTICAL  RROIAL 


66 


REFRRCT 

APPLIED 


ON  LINE  A 1  STN  =  3500  1 
BANDPASS  FILTER  5-0  — 


5.0  HZ 


5.45  TINE  IN  SEC.  5.95 


A 


IZ  .  ( R WRY  ) 


Figure  5.13  Plots  of  four  components  of  ground  motion  and 
radial  versus  vertical  ground  motion.  The  apparent  angle  of 
emergence  from  horizontal  boundary  is  31°. 
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REFRACTION  LINE  A2  STNr21005 
APPLIED  BANDPASS  FILTER  5-0—15.0  HZ 


12.00  TIME  IN  SEC.  12.50 
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Figure  5.14  Plots  of  four  components  of  ground  motion  and 
radial  versus  vertical  ground  motion.  The  apparent  angle  of 
emergence  from  horizontal  boundary  is  42c. 
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REFRACTION  LINE  N1  STN  =  35002 
APPLIED  BANDPASS  FILTER  5.0—15.0  HZ 


21  .20  TIME  IN  SEC .  21  .70 


Figure  5.15  Plots  of  four  components  of  ground  motion  and 
radial  versus  vertical  ground  motion.  The  apparent  angle  of 
emergence  from  horizontal  boundary  is  50°. 
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REFRRCT ION  LINE  F2  STN=21004 
APPLIED  BANDPASS  FILTER  5-0— 15-0  HZ 


17.00 


TIME  IN  SEC.  17.50 


Figure  5.16  Plots  of  four  components  of  ground  motion  and 
radial  versus  vertical  ground  motion.  The  apparent  angle  of 
emergence  from  horizontal  boundary  is  51°. 
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where : 

y  =  Poisson ' s  ratio, 

fi  =  angle  of  emergence  measured  from  the  horizontal. 

Cumming,  Clowes  and  Ellis  (1979)  concluded  that  Poisson's 
ratio  in  British  Columbia  was  0.23.  This  value  was  adopted 
to  compute  the  true  angle  of  emergence.  This  Poisson's  ratio 
may  be  in  error  when  it  is  applied  to  data  from  Vancouver 
Island.  Mereu  (1966)  showed  that  an  increase  of  Poisson's 
ratio  from  0.2  to  0.3  had  the  effect  of  increasing  the 
apparent  angle  of  emergence  a  few  degrees  (the  maximum  was 
about  7  degrees).  Thus  we  conclude  that  Poisson's  ratio  is 


not  a  critical  factor 

in  the 

determination  of  angle 

of 

emergence  in  our  case. 

By  knowing  the  velocity 

of 

the  wave 

at  the  bottom 

of 

the  trajectory  from 

the 

travel-time 

analysis,  we 

can 

estimate  the  near  surface  velocity  by  Snell's  law. 

sin(i)  =  V0/V, 

where : 

i  =  the  angle  of  emergence  with  respect  to  the  normal, 

V o  =  near  surface  velocity, 

V ,  =  velocity  at  the  maximum  depth  of  penetration. 

The  apparent  angles  of  emergence  from  all  particle  motion 
diagrams  were  measured,  and  then  converted  to  true  angles  of 
emergence  measured  from  the  vertical  in  order  to  compute  the 
near  surface  velocities  (Table  3). 

From  Table  3  the  surface  velocity  fluctuates  from  4.0 
to  5.27  km/s  in  this  particle  motion  analysis.  The 
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line 

geophone 
stat ion 

measured 

azimuth 

( degree ) 

angle  of 
emergence 

( degree ) 

velocity 
particle 
analy s i s 
km/s 

veloc i ty 
plane 
layer 
km/s 

A  1 

35001 

40 

56.3 

5.27 

5.28 

A2 

21005 

300 

43.4 

4.4 

4.57 

N  1 

35002 

307 

35.3 

4 . 0 

— 

F2 

21004 

1  32 

35.0 

4.0 

5.33 

Table  3....  Near  surface  velocities  from  the  particle  motion 
analyses  compared  with  those  from  plane-layered 
interpretation . 


uncertainty  of  Poisson's  ratio,  error  of  the  azimuth,  and 
the  assumption  of  a  homogeneous,  isotropic,  elastic  medium 
definitely  contribute  certain  errors  in  obtaining  the 
surface  velocities  at  various  locations.  When  the  errors  and 
assumption  are  taken  into  consideration,  the  surface 
velocities  from  particle  motion  analysis  are  considered  to 
agree  with  velocities  from  the  travel-time  analysis  based  on 
the  travel  time  to  the  first  detector.  Since  the  available 
data  are  limited  and  the  surface  velocities  from  the 
travel-time  analysis  are  not  well  defined,  it  may  be 
justified  to  use  the  velocity  of  4.0  km/s  from  particle 
motion  analysis  as  a  lower  bound  for  the  surface  velocity 
and  5.33  km/s  from  travel-time  data  as  a  upper  bound. 


6.  INTERPRETATION 


6 .  1  Introduction 

The  prime  objective  of  this  study  is  to  obtain  crustal 
structural  information  under  Vancouver  Island. 
Interpretations  of  crustal  structure  in  this  thesis  are 
based  on  the  analysis  of  refraction  data  through  the 
application  of  horizontal,  plane-layered  models,  and  through 
the  inversion  of  travel-time  curves  by  the  WHB  integral  and 
includes  synthetic  modelling  of  travel-time  branches.  The 
interpreted  results  of  four  crustal  seismic  sections 
evaluated  by  each  technique  are  combined  to  provide  two 
preliminary  models  in  the  study  area.  Comparison  of  the  two 
methods  indicates  that  the  model  obtained  by  the  WHB 
integral  technique  seems  to  represent  a  more  realistic  model 
than  the  horizontal-layered  interpretation. 


6.2  A  Preliminary  Model  from  Plane-layered  I nt repretat ion 

By  combining  the  results  of  the  individual  seismic 
sections  from  the  plane-layered  interpretation,  a  tentative 
regional  crustal  model  may  be  constructed  for  the  study 
area.  Because  of  limitations  of  the  data  and  the  large 
separation  distances  between  recording  sites,  this  model 
should  be  regarded  as  speculative  at  this  time. 

Fig. 6.1  shows  the  inferred  crustal  model.  This  consists 
of  4  horizontal  layers.  The  low  velocity  surface  layer  is 
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Figure  6.1  A  velocity  model  from  plane-layered 
interpretation  along  the  profile  lines  in  Fig. 3.1.  Contours 
are  velocity  in  km/s.  Shot  points  are  marked  with  arrow. 
Reliability  of  the  contours  is:  solid  line=good,  dashed 
1 ines=poor . 
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thin  towards  the  centre  of  the  profile,  and  become  thicker 
at  both  ends.  The  maximum  thickness  of  "sediments"  is  about 
2.1  km  at  the  south  end  and  the  minimum  is  0.51  km  at  the 
centre . 

The  depth  as  well  as  velocity  above  the  7.0  km/s 
refractor  shows  quite  substantial  changes  across  the 
profile.  The  velocity,  proceeding  from  both  ends  of  the 
profile  towards  the  centre,  decreases  from  approximately  6.5 
km/s  to  6.35  km/s,  and  there  is  a  significant  increase  of 
thickness  towards  the  centre.  These  lateral  variations  of 
velocities  and  thicknesses  may  suggest  a  disturbed  zone  or  a 
zone  of  major  changes  in  rock  compositions  across  the 
central  profile.  On  the  whole,  the  upper  crust  is  relatively 
homogeneous  except  at  the  centre  of  the  profile. 

The  observed  discontinuity  of  Pn  arrivals  certainly 
obscures  the  interpretation  of  the  Moho  arrivals.  This 
discontinuity  can  be  interpreted  as  indicating  the  existence 
of  a  low  velocity  zone  beneath  the  7.0  km/s  refractor.  If 
such  a  low  velocity  zone  is  present,  the  interpretation 
based  on  plane-layered  models  may  produce  serious  error  in 
estimating  the  Moho  depth,  because  it  violates  the  basic 
principle  of  plane-layered  modelling  in  which  velocity 
should  increase  with  depth.  However,  the  observations  of  Pn 
arrivals  on  profile  FI  provide  an  indirect  way  to  estimate 
the  depth  of  the  Moho.  Unfortunately,  the  high  noise  levels 
on  the  reversed  profile  make  the  picks  of  the  first  arrivals 
very  difficult,  thus  the  velocity  and  depth  of  the  Moho  are 
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not  well  constrained  by  the  refraction  data.  Since  we 
speculate  the  existence  of  the  low  velocity  zone,  the  Moho 
depth  (35  km)  by  plane-layered  interpretation  should  be 
regarded  as  tentative  at  this  time,  and  probably  represents 
a  minimum  crustal  thickness. 


6.3  A  Preliminary  Model  From  the  Synthetic  Seismogram  and 
WHB  Inversion  of  Travel  Times 

The  interpreted  results  from  the  WHB  inversion  of 
travel  times  and  synthetic  seismogram  modelling  were  also 
combined  to  form  a  tentative  interpretation  in  Fig. 6. 2.  The 
construction  of  this  model  was  based  on  the  following 
assumptions : 

1.  The  thicknesses  of  surface  "sediments"  were  taken  from 
the  plane-layered  interpretation. 

2.  The  velocity  increases  continuously  with  depth  to  about 
7.0  km/s  at  a  depth  of  approximately  17  km.  The  depth  to 
this  refractor  was  derived  by  combining  the 
velocity-depth  functions  of  four  seismic  sections. 

3.  The  WHB  inversion  of  travel-time  data  up  to  a  distance 
of  about  210  km  from  the  Fuca  profiles  yielded  a  maximum 
velocity  of  7.30  km/s  at  a  depth  of  27  km. 

4.  A  transition  zone  was  introduced  beneath  the  7.3  km/s 
refractor  by  synthetic  modelling.  The  velocity  in  the 
transition  zone  decreases  gradually  from  7.3  km/s  to  6.2 
km/s  at  the  depth  of  38  km. 
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Figure  6.2  A  velocity  model  from  synthetic  seismogram  and 
WHB  inversion  of  travel  times  along  the  profile  lines  in 
Fig. 3.1.  Contours  are  velocity  in  km/s.  Shot  points  are 
marked  with  arrow.  Reliability  of  the  contours  is:  solid 
line=good,  dashed  lines=poor. 
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5.  The  low  velocity  of  6.2  km/s  stays  constant  for  a 
thickness  of  6  km  and  sharply  jumps  to  a  velocity  of  8.0 
km.  This  low  velocity  channel  was  introduced  in  order  to 
account  for  the  large  amplitude  events  which  may  be  Moho 
reflections . 

This  model  in  Fig.  6.2  exhibits  a  different  crustal 
structure  compared  to  Fig.  6.1.  The  thickness  of  the  6.5 
km/s  refractor  is  fairly  constant  along  the  whole  profile  in 
Fig.  6.2.  There  are  slight  variations  of  thicknesses  south 
of  the  shot  point  A.  The  lack  of  structural  variation  may 
suggest  that  the  structure  down  to  the  6.5  km/s  refractor  is 
fairly  homogeneous. 

At  the  depth  about  17  km  the  7.0  km/s  refractor  is 
fairly  horizontal  at  both  ends  of  the  profile,  but  the  depth 
varies  substantially  near  the  south  of  location  A,  in  which 
the  variation  of  depths  is  well  represented  by  an  apparent 
syncline.  The  difference  in  depth  determined  from  southern 
profile  to  the  northern  profile  is  up  to  4  km.  This 
indicates  that  the  ray  path  that  travels  from  A  to  F  is 
considerably  different  from  the  ray  path  that  travels  from  A 
to  N.  One  possible  explanation  to  this  large  discrepancy  in 
depth  is  due  to  a  local,  anomalous  structure  south  of  shot 
point  A.  However,  if  this  anomalous  structure  exists,  the 
refaction  data  in  this  study  do  not  provide  enough 
information  to  determine  the  extent  of  this  zone,  therefore 
this  interpretation  is  to  some  extent  arbitrary  because  of 
significant  lateral  inhomogen iety . 
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Besides  the  uncertainty  in  the  depth  to  the  7.0  km/s 
refractor,  the  observation  of  a  discontinuity  in  travel 
times  contributes  another  uncertainty  to  estimate  the  depth 
of  the  Moho.  Langston  (1981)  concluded  the  existence  of  a 
low  velocity  zone  beneath  Vancourer  Island.  This  evidence 
provides  further  support  to  our  present  model.  Although  the 
thickness  of  the  low  velocity  zone  and  the  depth  to  the  Moho 
(44  km)  were  found  by  synthetic  modelling,  both  should  be 
considered  to  be  speculative  at  the  present  time.  Further 
seismic  surveys  are  required  to  establish  a  well-defined 
depth  of  the  Moho. 


6.4  Comparison  of  Two  Models 

Comparsion  of  the  models  in  Fig.  6.1  and  Fig.  6.2  shows 
that  the  depth  to  each  refracting  layer  from  the 
plane-layered  model  is  less  than  the  depth  from  the  WHB 
inversion  of  travel  times.  For  example,  the  maximum 
difference  in  depth  to  the  7.0  km/s  refractor  is  up  to  about 
7.3  km.  This  result  seems  to  be  quite  reasonable,  for  the 
WHB  integral  takes  into  consideration  velocity  increases 
with  depth  and  this  should  give  a  thicker  section  than  that 
obtained  by  assuming  a  constant  velocity.  In  the  real 
geological  situation,  the  assumption  of  velocity  increasing 
with  depth  is  presumably  more  realistic  than  the  assumption 
of  constant  velocities  in  a  small  number  of  layers.  The 
observed  discontinuity  of  travel  times  on  the  Fuca  profiles 
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certainly  favors  the  model  in  Fig.  6.2  which  accounts  for 
the  low  velocity  zone  indicated  by  the  travel-time 
discontinuity.  In  addition,  a  structural  interpretation 
based  on  an  iterative  modelling  of  the  reflected  and 
refracted  P-wave  travel  times  and  amplitudes  (McMechan  and 
Spence,  1982)  indicates  a  crustal  model  similar  to  that  in 


Fig.  6.2. 

The  depths  to 

both 

6.5 

km/s  and 

7.0 

km/s 

refractors 

are  consistent 

for 

both 

models,  but 

the 

model 

from  McMechan  and  Spence  has  a  very  complex  structure  near 
the  centre  of  the  profile.  Based  on  this  evidence,  the  model 
obtained  by  combining  the  synthetic  seismogram  and  the  WHB 
inversion  of  travel  times  is  preferred  as  the  final  model  in 
this  study. 


7.  CONCLUSIONS 


The  limited  extent  and  the  nature  of  the  data  in  the 
present  investigation  restrict  any  conclusions  to  only 
provisional  comments.  The  overall  paucity  of  observations  of 
Pn  arrivals  for  Vancouver  Island  also  limits  any  possible 
conclusions  about  the  regional  crustal  structure.  However, 
if  the  features  of  Fig.  6.2  represent  the  general  structure 
for  the  region,  several  conclusions  can  be  made: 

1.  The  surface  "sediments"  (low  velocity  materials)  are 
thinner  towards  the  centre  of  the  profile. 

2.  The  upper  crust  to  a  depth  of  about  5  km  is  well 
constrained  by  the  reversed  refraction  data  and  has  a 
homogeneous  velocity. 

3.  The  depth  to  the  7.0  km/s  refractor  shows  substantial 
change  of  crustal  thicknesses  to  the  south  of  location 
A.  This  suggests  the  existence  of  a  locally  anomalous 
structure  in  that  region. 

4.  Although  refracted  arrivals  from  the  Moho  have  been 
identified,  its  depth  and  shape  remain  uncertain, 
however,  a  depth  of  44  km  with  a  mantle  velocity  near 
8.0  km/s  is  reasonably  consistent  with  the  data. 

The  refraction  data  in  the  present  study  provide 
indirect  information  as  to  the  crustal  thickness  beneath 
Vancouver  Island.  In  attempting  to  interpret  the 
refraction  data  to  obtain  the  lower  crustal  and  upper 
mantle  structures  under  Vancouver  Island,  we  feel  that 
the  interpretation  of  the  present  data  is  somewhat 
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arbitrary  and  less  reliable  than  would  be  desirable.  The 
limitations  in  this  case  are  mainly  due  to  a  lack  of 
seismic  control  along  the  refraction  profiles 
particularly  at  large  distances,  and  the  complex 
structure  beneath  Vancouver  Island.  Larger  shots  and 
closer  spacing  of  recording  sites  are  suggested  to 
obtain  reliable  Moho  arrivals  for  further  studies.  If  a 
low  velocity  zone  exists  at  depth,  even  a  more  detailed 


refraction 

interpretation 

may 

be  ambiguous 

.  An 

indicat  ion 

of  complex  structure 

particularly 

in  the 

centre  of 

the  Island  and 

the  possibility  of 

a  low 

veloc i ty 

zone  at  depth 

might 

be  resolved 

by  a 

near-vertical  reflection  study. 
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APPENDIX  - COMPUTER  PROGRAMS 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THIS  PROGRAM  COMPUTES  THE  AMPLITUDE  AND  PHASE  SPECTRA 
OF  UNIVERSITY  OF  ALBERTA  RECORDING  SYSTEM  AND 
SEISMOMETER  IN  1980  SEISMIC  PROJECT  IN  VANCOUVER 
ISLAND. 

FQ  =  FREQUENCY  VECTOR  OF  THE  SPECTRUM 

AMP  =  AMPLITUDE  VECTOR  OF  THE  SPECTRUM 

PHASE  =  PHASE  ANGLE  VECTOR  OF  THE  SPECTRUM 
TRANS  =  TRANSFER  FUNCTION  OF  AMPLIFIER  AND 

SEISMOMETER  IN  S-PLANE  REPRESENTATION 


DIMENSION  FQ (81) , BLOGF ( 100) ,  AMP  ( 100) , BLOGA ( 100) 
DIMENSION  TITLEX( 60 ) ,TITLEY( 60 ) 

DIMENSION  BLOGAM( 1 00) ,PHASE( 1 00) 

COMPLEX  A, I ,S, TRANS 

COMMON  TI TLEX , TI TLEY , NTI TLX , NTI TLY , AXLEN , AYLEN 


SEE  THE  EXPLANATION  OF  VARIABLES  IN  SUBROUTINE  PLOTA . 

READ (5,111)  TITLEX 
READ ( 5 ,111)  TITLEY 

111  FORMAT ( 6 0A4) 

WRITE ( 6 , 1 1 1 ) TITLEY 
READ (5,112) AXLEN , AYLEN 

1 1 2  FORMAT ( 2F6 . 0 ) 


READ (5,113) NTI TLX , NTI TLY 
1 13  FORMAT (216) 

J.  -  A.  .  | ^  ■*.  tl .  ^  tl*  ^  ^  ^  *1*  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  .1*  ^  ^  ^  ^  ‘y* 

^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^ 

I=(0.0, 1 .0) 

ARC=6. 283 185 


FMI N=THE  MINIMUM  FREQUENCY 
FMAX=THE  MAXIMUM  FREQUENCY 
FACTOR=THE  GAIN  OF  THE  AMPLIFIER  (58DB) 

READ (5,4) FMI N , FMAX , FACTOR 
4  FORMAT (3F9. 4) 

NN=  1 

FREQ=FMI N 
100  S=I*FREQ 

FQ (NN ) =FREQ 

TRANS =  THE  TRANSFER  FUNCTION  OF  AMPLIFIER  +SEI SMOMETER 

A=TRANS ( S ) 

AMP ( NN ) =CABS ( A ) 

PHASE ( NN ) =ATAN2 ( AIMAG ( A ) , REAL (A) )* 360. /ARC 
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BLOGF ( NN ) =ALOG 1 0 ( FREQ ) 

CONVERT  AMPLITUDE  INTO  DECIBEL 

BLOGA ( NN ) = ALOG 1 0 ( AMP ( NN ) )*20 . O+FACTOR 
BLOGAM ( NN ) =  ALOG 1 0 ( AMP ( NN ) ) 

WRI TE ( 6 , 1 8 ) FREQ , AMP ( NN )  , BLOGAM ( NN ) , PHASE ( NN ) 

18  FORMAT (4( IX, FI  1 .4)  ) 

FREQ=FREQ* 1 . 1 547820 
NN=NN+ 1 

IF  ( FREQ . LE . FMAX ) GO  TO  100 
N1=NN- 1 

CALL  PLOTA (BLOGF, BLOGAM, N1 ) 

STOP 

END 

COMPLEX  FUNCTION  TRANS (S) 

COMPLEX  S 
TRANS= ( 1 ) 

1  ( (S/. 1 0 )/( 1 .0+S/. 10) )*1 .0/( 1 .0+1 . 2*S/4700 . 

THE  TRANSFORMER 

1  +( S/4700. )**2)*1 . 0/ ( 1 . O+S/53 . 1 ) 

THE  PREAMPLIFIER 
1  *(S/. 1 00 )/( 1 .0+S/. 100) 

THE  D.C. OFFSET  FILTER 
1  *  1 . 0/( 1 . O  +  S/53 .  1 ) 

THE  GAIN  CONTROL 

1  *105. /( 105. +105. *2. 1 105*S/50.+45.*(2. 1 105*S/50. )**2 
1  +10.* (2. 1 1 05*S/50 . ) *  *  3+ ( 2  .  1  1  05*S/50. )**4) 

THE  FOUR  POLE  BESSEL  FILTER  WITH  A  FUDGE  FACTOR 
FOR  FREQUENCY 

1  * ( S/2 . 0 ) *  *  3/ ( 1 ,0+.76*S/2.0+(S/2.0)**2)*8. 33*10.0/ 

1  ( 10. 0+4)*2*3. 1416*2 

THE  SEISMOMETER 
RETURN 
END 


SUBROUTINE  PLOTA (X,Y,N) 

REQUIRED  INPUTS 

THE  TITLES  OF  X-  AND  Y-  AXIS 

AXLEN-AYLEN:  DESIRED  LENGTH  OF  X  AND  Y  AXIS  IN  INCHES 
NTITLX-NTITLY:  NUMBER  OF  CHARS  IN  X-AND  Y-TITLES 
X-Y:  THE  VECTORS  TO  BE  PLOTTED 
NDATA :  NUMBER  OF  ELEMENTS  IN  VECTOR  X  OR  Y 

DIMENSION  X ( N ) , Y ( N ) 

DIMENSION  TITLEX(60) ,TITLEY(60) 

COMMON  TI TLEX , TI TLEY , NTI TLX , NTI TLY , AXLEN , AYLEN 

PLOT  EVERY  DATA  POINT 
I REP= 1 


CALL  PLOTS 
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ENLARGE  OR  REDUCE  THE  SIZE  OF  THE  ENTIRE  PLOT 
CALL  FACTOR (0.8) 

MOVE  ORIGIN  IN  A  BIT  TO  ALLOW  FOR  LABELS 
CALL  PLOT ( 1 .0,3. 0,-3) 

ENLARGE  OR  REDUCE  THE  SIZE  OF  THE  ENTIRE  PLOT 

SCALE  THE  INPUT  VECTORS 

CALL  SCALE ( X , AXLEN , N , I  REP ) 

CALL  SCALE ( Y , AYLEN , N , I  REP ) 


X ( N+ 1 )=- 1 . 0 
X ( N+2 ) =0 . 5 
Y(N+1 ) =0 . 0 
Y ( N+2 ) = 1 .3333 

PLOT  X-AND  Y-  AXIS 

CALL  LOGAX ( 0 . 0 , 0 . 0 , TI TLEX , -NT I  TLX , AXLEN , 0 . 0 , -  1 . 0 , 
&  -3. 0,3.0) 

CALL  LOGAX ( 0 . 0 , 0 . 0 , TI TLEY , NTI TLY , AYLEN, 90. ,-0.0, 

&  -4. 0,3.0) 

CALL  LINE ( X , Y , N , 1,0,3) 

TERMINATE  THE  PLOT 

CALL  PLOT( 0 . 0 , 0 . 0 , 999) 

RETURN 

END 
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C 

C 

C  THIS  PROGRAM  COMPUTES  THE  DEPTHS  OF  N  LAYERS  WITH 
C  GIVEN  INPUTS  : 

C  VELOCITIES  AND  INTERCEPT  TIMES  OF  GIVEN  LAYERS. 

C 

C  NL  =  NUMBER  OF  LAYERS  TO  BE  COMPUTED  NOT  INCLUDE 
C  HALF  SPACE 

C  NV  =  NUMBER  OF  VELOCITIES 

C  TAU  =  THE  INTERCEPT  TIMES 

C  AW  =  VECTOR  TO  STORE  THE  TITLE  OF  THE  LINE 

C  V  VELOCITIES  OF  N  LAYERS 

C  H  THICKNESS  OF  INDIVIDUAL  LAYER 

C 
C 

REAL  TAU( 10) ,V( 10) ,H( 10) ,AW(30) 

DATA  TAU/ 10*0. 0/,V/ 10*0. 0/,H/ 10*0.0/ 

READ (5, 4) (AW (I ) ,1=1 ,30) 

WRI TE (6,4) (AW(I )  ,1  =  1  ,30) 

4  FORMAT ( 30A4 ) 

READ  (5,2)  NV 
2  FORMAT (II) 

NOTE  TAU ( 1 )  MUST  BE  0.0 

READ ( 5 ,  1 0 ) ( TAU ( I  )  ,  I  =  1  , NV ) 

10  FORMAT ( 7F 1 0 . 3 ) 

READ (5,20)(V(I),I=1,NV) 

20  FORMAT ( 8F9 . 3 ) 

DO  90  N  =  2 , NV 
HSUM=0 . 0 
NLAYER=N-2 

I F ( NLAYER . LE . 0 )  GO  TO  110 
DO  100  K= 1 , NLAYER 
VI =SQRT (V(N)**2~V(K)**2) 

V2=SQRT (V(N)**2-V(N~1 )**2) 

HSUM=HSUM+H ( K ) * ( V ( N~ 1 )*V1 )/(V(K)*V2) 

WRI TE ( 6 , 1 0 1 ) N , K , V 1 , V2 , HSUM 
101  FORMAT ( 21 5, 4F1 0.3) 

100  CONTINUE 
C 

110  V3=SQRT(V(N) **2-V(N- 1 )**2) 

C 

H(N-1 )=TAU(N)*V(N-1 ) *V(N)/( 2 . 0*V3 ) -HSUM 
90  CONTINUE 
C 

WRITE (6, 30) 

3  0  FORMAT ( ' 0 '  , 4X, ' INTERPT  TI ME '  , 5X , ' VELOCI TY '  , 5X  , 

&  ’THICKNESS' ,7X, ’DEPTH’ ) 

DEPTH=0 . 0 
DO  40  L= 1 , NV 
DEPTH=DEPTH+H(L) 

IF(L.EQ.NV)  DEPTH=0 . 0 
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WRITE (6, 45)  TAU(L) ,V(L) ,H(L) , DEPTH 
45  FORMAT ( ’0',5X,F10.3,5X,F8.3,5X,F9.3,5X,F9.3) 
40  CONTINUE 
STOP 
END 
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THIS  PROGRAM  COMPUTES  THE  VELOCITY-DEPTH 
FUNCTION  FROM  THE  INVERSION  OF  TRAVEL  TIMES 
BY  WHB  INTEGRAL,  USING  A  LOW  ORDER  POLYNOMIAL 
TO  FIT  WITH  THE  OBSERVED  TRAVEL-TIME  DATA 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


8 


1 1 
1  3 
14 
1  0 


25 

20 


DIMENSION  A(20) , B ( 20 ) , DI ST ( 100) , TIME ( 100) ,X( 100) 
DIMENSION  C( 10) ,S(20) , TY (100), AW (30) ,VEL( 100) 
DIMENSION  Z (  50 0 )  ,VMAX(  50  0) , Z 1  ( 500 ) , Y ( 100) 
DIMENSION  RX ( 3 ) , RT ( 3 ) 

COMMON  AW 

REAL*8  P ( 200 ) ,T(200) 

DATA  C/1  0*0.0/ 


AW  =  VECTOR  TO  STORE  THE  TITLE  OF  REFRACTION  LINE 
RX  =  DISTANCES  USED  TO  REDUCE  DATA  TO  THE  BASE 
OF  THE  SURFACE  LAYER 

RT  =  TIMES  USED  TO  REDUCE  DATA  TO  THE  BASE  OF  THE 
SURFACE  LAYER 

DIST  =  DISTANCES  OF  TRAVEL-TIME  DATA 
TIME  =  TIMES  OF  TRAVEL-TIME  DATA 

TY  =  TIMES  COMPUTED  FROM  POLYNOMIAL  COEFFICIENTS 
Z  =  DEPTHS  CALCULATED  FROM  WHB  INTEGRAL 
VEL  =  VELOCITIES  CALCULATED  FROM  WHB  INTEGRAL 
N  =  INPUT  NUMBER  OF  OBSERVATIONS  OF  X  AND  Y 
MD  =  MAX  DESIRED  DEGREE  OF  THE  POLYNOMIAL  FIT 
ID  =  OUTPUT  DEGREE  OF  THE  FITTED  MODEL 
C  =  OUTPUT  COEFFICIENTS  OF  POLYNONIMAL  FIT  IN 
ASCENDING  ORDER 

P  AND  T  =  WORKING  SPACE  MUST  BE  DOUBLE  PRECISION 


RSQ= 100. 

READ (5, 8) (AW (I), 1=1, 30) 

WRI TE ( 6 , 8 ) ( AW ( I ) , I  =  1 , 3  0 ) 

FORMAT ( 3 0A4 ) 

READ  (5, 1 1 ) (RX(K) ,K=1 , 3) 

READ  (5, 1 1 ) (RT(K) ,K=1 ,3) 

FORMAT (3F1 0.3) 

WRI TE ( 6 , 13) (RX(K) ,K=1 ,3) 

FORMAT ( 'O'  ,  '  RX' , 3F10.3) 

WRI TE ( 6 , 14) (RT(K) ,K=1 , 3) 

FORMAT ( ' O' , '  RT' , 3F10.3) 

READ (5,10)  MD , N 
FORMAT ( 21 6 ) 

DO  20  1=1 ,N 

READ (5,25)  DI ST ( I ) , TIME ( I ) 
FORMAT ( 2F 12.4) 

CONTINUE 

CALL  REDUTI (N, DIST, TIME, RT,RX) 
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DO  30  1= 1 ,N 
X(I ) =DI ST ( I ) 
30  Y ( I ) =TIME ( I ) 
TIME ( 1 ) =0 . 0 


RLFOTH  FROM  IMSLLIB  —  FIT  A  UNIVARIATE  CURVILINEAR 
REGRESSION  MODEL  USING  ORTHOGONAL  POLYNOMIALS 

CALL  RLFOTH ( DI ST , Y , N , RSQ , MD , I D , P , C , S , A , B , I ER ) 

RLDOPM  FROM  IMSLLIB  —  TRANSFORMS  A  FITTED  ORTHOGONAL 
POLYNOMIAL  MODEL  INTO  A  POLYNOMIAL  FUNCTION  OF  THE 
ORIGINAL  INDEPENDENT  VARIABLE 

CALL  RLDOPM (C,ID,A,B,T) 

THE  OUTPUT  OF  C  CONTAINS  THE  CONSTANT  TERM 
AND  COEFFICIENTS 

LD= I D+ 1 
DO  40  K=  1  ,  N 
CONST=C ( 1 ) 

YY=CONST 
DO  50  KK=  2 , LD 

50  YY=YY+C (KK) *X(K) ** (KK~ 1 ) 

TY ( K ) =YY 
40  CONTINUE 

WRI TE ( 6 , 60 )  ID 

60  FORMAT ( ' 0 T  , 5X,  ' DEGREE  OF  POLYNOMIAL  FIT  I  S’, I  3) 

WRI TE (6,65) 

65  FORMAT ( ’ 0 ’ , 1 5X, ' THE  COEFFS  OF  POLY.  FITTED  MODELS') 
WRI TE ( 6 , 7  0 ) ( C ( K ) , K= 1 , LD ) 

70  FORMAT ('O' , 5 ( 3X , E 1 4 . 6 ) ) 

WRITE (6, 80) 

8  0  FORMAT ( ' 0 '  , 5X ,  ' DI STANCE '  ,  1  OX , ' TI ME '  , 5X , 

&  'THEORETICAL  TI ME '  ,  1  OX ,  ' DI FF '  ) 

STD=0 . 0 

DO  90  K= 1 ,N 

DIFF=TY(K)-Y(K) 

STD=STD+DIFF**2 

WRI TE ( 6 , 95 )  X(K) ,Y(K) ,TY(K) , DI FF 
95  FORMAT ('O'  , 5X , F9 . 3 , 5X , F9 . 3 , 1  OX , F9 . 3 , 5X , F9 . 3 ) 

90  CONTINUE 

STD=STD/N 
WRITE (6, 91)  STD 

91  FORMAT ( ’ 0 ', 4X, ' SQUARE  NORM  OF  TIME  IS',F9.6) 

CALL  PLOTS 

CALL  PLOTA ( X , TIME , TY , N ) 
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C 

N  =  4  3 

DO  93  K=  1  ,  N 
93  X ( K ) = ( K- 1 ) *5. 0 
C 

CALL  Z VEL ( X  ,  N  ,  0 . 2  ,  C  , LD , Z , VMAX ) 

WRI TE ( 6 , 154) 

1 54  FORMAT ( 8X , ' DI STANCE ’ , 8X , ’ DEPTH ' , 8X , ? VMAX ’ ) 

DO  150  K= 1 , N 

WRI TE ( 6 , 155)  X ( K ) ,  Z  ( K ) , VMAX ( K ) 

155  FORMAT ('O' , 4X , 3F 1 2 . 4 ) 

150  CONTINUE 

CALL  PLTA ( Z , VMAX , N ) 

C 

STOP 

END 

C 

SUBROUTINE  PLOTA ( X , Y , TY , N ) 

C 

C  THIS  SUBROUTINE  PLOTS  THE  TIME-DISTANCE  CURVE. 

C 

C  REQUIRED  INPUTS 

C  THE  TITLES  OF  X-  AND  Y-  AXIS 

C  AXLEN-AYLEN:  DESIRED  LENGTH  OF  X  AND  Y  AXIS 

C  NTITLX-NTITLY:  NUMBER  OF  CHARS  IN  X-AND  Y-TITLES 

C  X-Y:  THE  VECTORS  TO  BE  PLOTTED 

C  NDATA :  NUMBER  OF  ELEMENTS  IN  VECTOR  X  OR  Y 

C 

DIMENSION  X ( 2 ) ,  Y  (  2  ) , TY ( 2 ) 

COMMON  AW 

PLOT  EVERY  DATA  POINT 
I REP  = 1 
AXLEN=5 . 0 
AYLEN=4 . 0 

ENLARGE  OR  REDUCE  THE  SIZE  OF  THE  ENTIRE  PLOT 
CALL  FACTOR ( 1.0) 

MOVE  ORIGIN  IN  A  BIT  TO  ALLOW  FOR  LABELS 
CALL  PLOT ( 1.0,  11. 0,-3) 

SCALE  THE  INPUT  VECTORS 

X(N+1 ) =0 . 0 
X ( N+2 ) =44 . 0 
Y(N+1 ) =0 . 0 
Y (N+2 ) = 1 0 . 0 

PLOT  X-AND  Y-  AXIS 

CALL  AXIS (0.0, 0.0,  'DISTANCE  KM'  , -  1 1 , AXLEN , 0 . 0 , 

&  X(N+1 ) , X ( N+2 ) ) 

CALL  AXIS(0. 0,0.0, 'TIME  SEC ’ , 8 , AYLEN , 90 . , Y ( N+ 1 ) , 
&  Y (N+2 ) ) 
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CALL  SYMBOL (0.4, 4. 4, 0.16, AW, 0.0, 30) 

C 

CALL  LINE ( X  ,  Y  ,  N  ,  1,0,3) 

C  CALL  LI NE ( X , TY , N ,  1  ,0,3) 

WRI TE ( 6 , 10)  X(N+1 ) ,X(N+2) ,Y(N+1 ) ,Y(N+2) 

10  FORMAT (’  SCALE  FACTOR  '  , 4F  15.2) 

TERMINATE  THE  PLOT 

CALL  PLOT (0.0, -4. 6, -3) 

RETURN 
END 

SUBROUTINE  PLTA(X,Y,N) 

THIS  SUBROUTINE  PLOTS  THE  VELOCITY-DEPTH  CURVE 


DIMENSION  X ( 2 ) , Y ( 2 ) 

PLOT  EVERY  DATA  POINT 
I REP= 1 


SCALE  THE  INPUT  VECTORS 
AXLEN=  5 . 0 
AYLEN=  3 . 0 
X ( N+ 1 ) =0 . 0 
X (N+2 ) =6 . 0 
Y(N+1 ) =5 . 0 
Y(N+2)  =  1  .0 

PLOT  X-AND  Y-  AXIS 

CALL  AXIS (0.0, 0.0, 'DEPTH  KM ' , -8 , AXLEN , 0 . 0 , 

&  X ( N+ 1 ) ,X(N+2 ) ) 

CALL  AXIS (0.0, 0.0, 'VELOCITY  KM/S ’ , 1 3 , AYLEN , 90 . , 
&  Y(N+ 1 ) , Y(N+2 ) ) 

CALL  L I NE (X,Y,N, 1 ,0,3) 

TERMINATE  THE  PLOT 

CALL  PLOT(0. 0,0. 0,999) 

RETURN 

END 


SUBROUTI NE  ZVEL ( X , LX , DELX , C , LC , Z , VEL ) 

THIS  SUBROUTINE  EVALUATES  THE  WHB  INTEGRAL  AT  ANY 
ARBITRARY  DISTANCE  STARTING  WITH  A  INITIAL  VELOCITY 
OF  5.30  KM/S. 

DIMENSION  X(2) ,C(2) ,Z(2) , VEL ( 2 ) 

ARCOSH ( Y ) =ALOG ( Y  +  SQRT ( Y*  *  2- 1 . 0) ) 
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K 1  =  1 

DO  1 0  KK= 1 , LX 
ORI G=X ( 1 ) 

XMAX=X ( K 1 ) -ORI G 
IF(XMAX.LE.O.O)  GO  TO  15 
SLOPE=DERV ( C , LC , XMAX ) 

VEL ( K 1 )= 1 /SLOPE 
V=VEL ( K 1 ) 

N=XMAX/DELX 
ZZ  =  0 . 0 

X 1 =xmax/n 

XX  =  X1 

DO  20  K= 1 , N 
SLOPE=DERV ( C , LC , XX) 
ARG=V*SLOPE 
IF  (K.EQ.N)  ARG=1.0 
ZZ=ZZ+ARCOSH(ARG) *DELX 
20  XX=XX+X1 
GO  TO  30 

15  VEL (K1) =5. 30 
Z ( K 1 ) =0 . 0 
K  1  =K 1  +  1 
GO  TO  10 

30  Z(K1 ) = 1 / 3 . 1416 *ZZ 
K  1  =K 1  +  1 
10  CONTINUE 
RETURN 
END 


REAL  FUNCTION  DERV(C,LC,X) 

THIS  FUNCTION  CALCULATES  THE  DERIVATIVE  AT  ANY 
POINT  FROM  THE  POLYNOMIAL  FUNCTION. 

DIMENSION  C ( LC ) 

Y=C ( 2 ) 

DO  10  L=  3 , LC 

10  Y=Y+C(L)*(L-1 )*X**(L-2) 

DERV=Y 

RETURN 

END 


SUBROUTI NE  REDUTI ( LX , X , T , RT , RX ) 

THIS  SUBROUTINE  REDUCES  THE  TRAVEL-TIME  DATA 
DOWN  TO  THE  BASE  OF  THE  SURFACE  LAYER 

DIMENSION  X ( 2 ) ,  T  (  2  ) , RT ( 2 ) , RX ( 2 ) , V ( 2 ) ,DEPTH(2) 
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IF  ( RT ( 3 ) . LE .0.0)  GO  TO  100 
R 1 = ( RT ( 1 ) -RT ( 2 ) ) / ( RX ( 2 ) -RX ( 1 ) ) 
R2=(RT(3)-RT(2) ) / ( RX ( 3 ) -RX ( 2 ) ) 

C 

T  (  1  )  =0 . 0 
DO  10  1=2, LX 
TORI G  =  T ( I  ) 

I F ( X ( I ) . LE . RX ( 2 ) )  T 1 =R 1  * ( RX ( 2 ) -X ( I ) ) +RT ( 2 ) 
I F ( X ( I ) . GT . RX ( 2 ) )  T1=R2*(X(I )-RX(2) )+RT(2) 
T ( I ) =T ( I )-T1 

WRITE (6, 20)  X(I ) ,TORIG,T(l ) ,T1 
20  FORMAT ('0' ,  '  H  1  ’ ,4(F10.4,5X)) 

10  CONTINUE 
GO  TO  110 
C 

100  R 1 = ( RT ( 2 ) -RT ( 1 ) ) / ( RX ( 2 ) -RX ( 1 ) ) 

T ( 1 ) =0 . 0 
C 

DO  40  I =2, LX 

I F ( X ( I )  . LE . RX ( 2 ) )  T 1 =R 1 *X ( I  )+RT( 1 ) 

T ( I ) =T ( I )-T1 

C  WRITE (6, 45)  X(l),T(l),T1 

C  45  FORMAT ( ’ O’ , ’HI’ , 3F12.4) 

40  CONTINUE 
1 1 0  RETURN 
END 


